Setup

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
# Set options for analysis
shutdown = FALSE

options(dplyr.print_max = 100, 
        brms.backend = 'cmdstan',
        brms.file_refit = 'on_change' # "on_change" or "never"
        )
# functions for kable printing. 
print_df <- function(df, width = "auto") {
  kable(df) %>%
    kable_styling(
      bootstrap_options = c("striped", "hover", "responsive"),
      full_width = F,
      position = "left",
      fixed_thead = T
    ) %>%
    column_spec(2:ncol(df), width = width) %>%
    row_spec(1:nrow(df), extra_css = "white-space: nowrap;") %>%
    scroll_box(width = '100%', height = '500px')
}

packing_ind <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 2, 6) %>%
  pack_rows(
    "Between-Person Effects", 7, 10) %>%
  pack_rows(
    "Random Effects", 11, 13) %>%
  pack_rows("Additional Parameters", 14, 18)
}

packing_apim <- function(df) {
  pack_rows(df,
    "Within-Person Effects", 3, 12) %>%
  pack_rows(
    "Between-Person Effects", 13, 20) %>%
  pack_rows(
    "Random Effects", 21, 26) %>%
  pack_rows("Additional Parameters", 27, 31)
}


# Function to prepare for export
prepare_for_export <- function(model) {
  rep <- model
  rep$coefs <- rownames(rep)
  final <- cbind(rep$coefs, model)
  return(final)
}

Data Preparation

Scales and Scores

predictors_to_center <- c("persuasion",
                          'pressure',
                          'pushing', 
                          'weartime',
                          'barriers',
                          'support',
                          'comf',
                          'reas')


predictors_not_to_center <- c('day', 'userID', 'coupleID', # independent variables that don't need to be centered.
                              'aff', # outcomes (not centered)
                              'pa_sub',
                              'pa_obj', 
                              'anticipAff', 
                              'reactance', 
                              'plan',
                              'isWeekend',
                              'studyGroup',
                              'got_JITAI',
                              'skilled_support'
                              )



all_variables <- c(predictors_not_to_center, predictors_to_center)


df <- df %>% 
  mutate(ss_affect1 = ss_affect1 + 1, # re-code items from 0-5 to 1-6 scale.
         ss_affect2 = ss_affect2 + 1, 
         ss_affect3 = ss_affect3 + 1, 
         ss_affect4 = ss_affect4 + 1, 
         
         aff = (ss_affect1 + ss_affect3 + (7 - ss_affect2) + (7 - ss_affect4)) / 4, #Affect Scale

         
         pa_sub = ifelse(ss_pa == 0, 0, ss_pa_min_solo + ss_pa_min_collab),
         
         pa_obj = ifelse(wear_minutes > 600, minutes_mvpa_non_filtered, NA),
         
         day = day/54,
         
         barriers = (ss_barr_1 + ss_barr_2 + ss_barr_3 + ss_barr_4 + ss_barr_5 + ss_barr_6 + ss_barr_7) / 7,
      
         plan = ifelse(
           is.na(ss_pa_no), 
           ifelse(ss_pa_yes_0 == 1 | ss_pa_yes_1 == 1 | ss_pa_yes_2 == 1 | ss_pa_yes_3 == 1, 1,0),
           ifelse(ss_pa_no ==1, 1, 0)
         ), 
         
         got_JITAI = ifelse(((userID %% 2 != 0) & (trig_target_slot_0 %in% c("A", "C") | trig_target_slot_1 %in% c("A", "C") | trig_target_slot_2 %in% c("A", "C"))) |
                              ((userID %% 2 == 0) & (trig_target_slot_0 %in% c("B", "C") | trig_target_slot_1 %in% c("B", "C") | trig_target_slot_2 %in% c("B", "C"))), 1, 0),
         
         got_JITAI = ifelse(is.na(got_JITAI), 0, got_JITAI),
         
         skilled_support = ifelse(
           (studyGroup == 3 & ((day * 54 + 1) > 28)) | (studyGroup != 3 & ((day * 54 + 1) > 7)), 
           1, 0
           )
         
         )%>% 
  rename(persuasion = sp_psc_more,
         pressure = sp_nsc_more,
         pushing = sp_push_plan,
         anticipAff = ss_affect_con,
         reactance = ss_reactance,
         weartime = wear_minutes,
         # aff_reactance_persuasion = sp_persuasion_aff,
         #aff_reactance_pressure = sp_nsc_aff, # explorative
         #aff_reactance_pushing = sp_push_plan_aff, # explorative
         
         support = sp_emo_pleasure,
         comf = sp_emo_comf,
         reas = sp_emo_reass
         ) 
# pushing can only occur if there was a plan! Therefore, we set it to zero if there was no plan. 
df$pushing[df$plan == 0] <- 0

center Data

prepare_df <- function(df) {
  df_wide <- df %>%
    # center
    bmlm::isolate(by = "userID",
                  value = predictors_to_center,
                  which = "both") %>%
    # distinguish
    group_by(userID) %>%
    mutate(person_mean_pa_obj = mean(pa_obj, na.rm = TRUE)) %>%
    ungroup() %>%
    
    group_by(coupleID) %>%
    mutate(
      role = case_when(
        person_mean_pa_obj == max(person_mean_pa_obj, na.rm = TRUE) ~ 'A', 
        TRUE ~ 'B'
    )) %>%
    ungroup() %>%
  
    # Reshape the data to wide format
    pivot_wider(
      id_cols = c(coupleID, day),
      names_from = c(role),
      values_from = -c(1, 2, 4)
    )
    
    
  # Double Stacking Outcomes 
    
  df1 <- df_wide
  df2 <- df_wide
  
  df1$is_A <- 1
  df1$is_B <- 0
  df1$AorB <- 1
  df1$pa_sub <- df1$pa_sub_A
  df1$pa_obj <- df1$pa_obj_A
  df1$aff <- df1$aff_A
  df1$reactance <- df1$reactance_A
  df1$anticipAff <- df1$anticipAff_A
  #df1$aff_reactance_pressure <- df1$aff_reactance_pressure_A
  #df1$aff_reactance_pushing <- df1$aff_reactance_pushing_A
  
  df2$is_A <- 0
  df2$is_B <- 1
  df2$AorB <- 2
  df2$pa_sub <- df2$pa_sub_B
  df2$pa_obj <- df2$pa_obj_B
  df2$aff <- df2$aff_B
  df2$reactance <- df2$reactance_B
  df2$anticipAff <- df2$anticipAff_B
  #df2$aff_reactance_pressure <- df2$aff_reactance_pressure_B
  #df2$aff_reactance_pushing <- df2$aff_reactance_pushing_B
  
  
  df_double <- rbind(df1, df2)

  
  # Double Stacking predictors
  # persuasion
  df_double$persuasion_self <- df_double$persuasion_A * df_double$is_A + df_double$persuasion_B * df_double$is_B
  df_double$persuasion_partner <- df_double$persuasion_A * df_double$is_B + df_double$persuasion_B * df_double$is_A
  
  df_double$persuasion_self_cw <- df_double$persuasion_cw_A * df_double$is_A + df_double$persuasion_cw_B * df_double$is_B
  df_double$persuasion_partner_cw <- df_double$persuasion_cw_A * df_double$is_B + df_double$persuasion_cw_B * df_double$is_A
  
  df_double$persuasion_self_cb <- df_double$persuasion_cb_A * df_double$is_A + df_double$persuasion_cb_B * df_double$is_B
  df_double$persuasion_partner_cb <- df_double$persuasion_cb_A * df_double$is_B + df_double$persuasion_cb_B * df_double$is_A
  
  #pressure
  df_double$pressure_self <- df_double$pressure_A * df_double$is_A + df_double$pressure_B * df_double$is_B
  df_double$pressure_partner <- df_double$pressure_A * df_double$is_B + df_double$pressure_B * df_double$is_A
  
  df_double$pressure_self_cw <- df_double$pressure_cw_A * df_double$is_A + df_double$pressure_cw_B * df_double$is_B
  df_double$pressure_partner_cw <- df_double$pressure_cw_A * df_double$is_B + df_double$pressure_cw_B * df_double$is_A
  
  df_double$pressure_self_cb <- df_double$pressure_cb_A * df_double$is_A + df_double$pressure_cb_B * df_double$is_B
  df_double$pressure_partner_cb <- df_double$pressure_cb_A * df_double$is_B + df_double$pressure_cb_B * df_double$is_A
  
  # pushing
  df_double$pushing_self <- df_double$pushing_A * df_double$is_A + df_double$pushing_B * df_double$is_B
  df_double$pushing_partner <- df_double$pushing_A * df_double$is_B + df_double$pushing_B * df_double$is_A
  
  df_double$pushing_self_cw <- df_double$pushing_cw_A * df_double$is_A + df_double$pushing_cw_B * df_double$is_B
  df_double$pushing_partner_cw <- df_double$pushing_cw_A * df_double$is_B + df_double$pushing_cw_B * df_double$is_A
  
  df_double$pushing_self_cb <- df_double$pushing_cb_A * df_double$is_A + df_double$pushing_cb_B * df_double$is_B
  df_double$pushing_partner_cb <- df_double$pushing_cb_A * df_double$is_B + df_double$pushing_cb_B * df_double$is_A
  
  # Weartime
  
  df_double$weartime_self <- df_double$weartime_A * df_double$is_A + df_double$weartime_B * df_double$is_B
  df_double$weartime_partner <- df_double$weartime_A * df_double$is_B + df_double$weartime_B * df_double$is_A
  
  df_double$weartime_self_cw <- df_double$weartime_cw_A * df_double$is_A + df_double$weartime_cw_B * df_double$is_B
  df_double$weartime_partner_cw <- df_double$weartime_cw_A * df_double$is_B + df_double$weartime_cw_B * df_double$is_A
  
  df_double$weartime_self_cb <- df_double$weartime_cb_A * df_double$is_A + df_double$weartime_cb_B * df_double$is_B
  df_double$weartime_partner_cb <- df_double$weartime_cb_A * df_double$is_B + df_double$weartime_cb_B * df_double$is_A
  
  # barriersiers
  df_double$barriers_self <- df_double$barriers_A * df_double$is_A + df_double$barriers_B * df_double$is_B
  df_double$barriers_partner <- df_double$barriers_A * df_double$is_B + df_double$barriers_B * df_double$is_A
  
  df_double$barriers_self_cw <- df_double$barriers_cw_A * df_double$is_A + df_double$barriers_cw_B * df_double$is_B
  df_double$barriers_partner_cw <- df_double$barriers_cw_A * df_double$is_B + df_double$barriers_cw_B * df_double$is_A
  
  df_double$barriers_self_cb <- df_double$barriers_cb_A * df_double$is_A + df_double$barriers_cb_B * df_double$is_B
  df_double$barriers_partner_cb <- df_double$barriers_cb_A * df_double$is_B + df_double$barriers_cb_B * df_double$is_A
  
  
  
  ## For Sensitivity Analyses necessary
  
  # Plans
  df_double$plan_self <- as.numeric(df_double$plan_A) * df_double$is_A + as.numeric(df_double$plan_B) * df_double$is_B
  df_double$plan_partner <- as.numeric(df_double$plan_A) * df_double$is_B + as.numeric(df_double$plan_B) * df_double$is_A
  
  # Support
  df_double$support_self <- df_double$support_A * df_double$is_A + df_double$support_B * df_double$is_B
  df_double$support_partner <- df_double$support_A * df_double$is_B + df_double$support_B * df_double$is_A
  
  df_double$support_self_cw <- df_double$support_cw_A * df_double$is_A + df_double$support_cw_B * df_double$is_B
  df_double$support_partner_cw <- df_double$support_cw_A * df_double$is_B + df_double$support_cw_B * df_double$is_A
  
  df_double$support_self_cb <- df_double$support_cb_A * df_double$is_A + df_double$support_cb_B * df_double$is_B
  df_double$support_partner_cb <- df_double$support_cb_A * df_double$is_B + df_double$support_cb_B * df_double$is_A
  
  # JITAIs
  df_double$got_JITAI_self <- df_double$got_JITAI_A * df_double$is_A + df_double$got_JITAI_B * df_double$is_B
  df_double$got_JITAI_partner <- df_double$got_JITAI_A * df_double$is_B + df_double$got_JITAI_B * df_double$is_A  
  
  
  df_double$pa_obj_log <- log(df_double$pa_obj)
      
  return(list(df_wide, df_double)) 
}


# And prepare original dataset
df_list <- prepare_df(df)

df_wide <- df_list[[1]]
df_double <- df_list[[2]]

Creating function to report measures

report_measures <- function(data, measures, ICC = TRUE) {
  
  if (ICC & any(sapply(data[,measures], is.factor))) {
    warning("Cannot add ICCs when Factors are in Dataframe.")
    ICC <- FALSE
  }
  
  measures_table <- as.data.frame(report::report_table(data[, measures]))

  measures_table <- measures_table %>%
    mutate(across(.cols = c('Mean', 'SD'),
                  .fns = ~format(round(.x, 2), nsmall = 2)))
  
  
  measures_table$Missing <- ifelse(
    is.na(measures_table$percentage_Missing),
    NA,
    paste0(round(measures_table$percentage_Missing), '%')
  )
  
  measures_table$Range <- ifelse(
    is.na(measures_table$Min) | is.na(measures_table$Max),
    NA,
    paste0(format(round(measures_table$Min, 2), nsmall = 2), '-', format(round(measures_table$Max, 2), nsmall = 2))
  )
    
  
  
  finished_df <- measures_table
  
  if (ICC) {
    cors <- wbCorr(data[,measures], data$coupleID)
    
    finished_df <- cbind(
      measures_table[, c('Variable', 'n_Obs', 'Missing', 'Mean', 'SD', 'Range')],
      ICC = format(round(get_icc(cors)$ICC, 2), nsmall = 2)
    )
  } else {
    measures_table$percentage_Obs <- ifelse(
      is.na(measures_table$percentage_Obs),
      NA,
      paste0(round(measures_table$percentage_Obs), '%')
    )
    
    finished_df <- measures_table[, c('Variable', 'Level', 'n_Obs', 'percentage_Obs', 'Missing', 'Mean', 'SD', 'Range')]
  }
  
  # Make them all Text for easier copying to excel. 
  finished_df[] <- lapply(finished_df, as.character)
  
  return(finished_df)
}

Sample Description

# mean of both relationship durations reported
df_wide$rel_duration_A <- df_wide$pre_rel_duration_y_A + df_wide$pre_rel_duration_m_A / 12 # in months
df_wide$rel_duration_B <- df_wide$pre_rel_duration_y_B + df_wide$pre_rel_duration_m_B / 12
df_wide$Relationship_Duration <- (df_wide$rel_duration_A + df_wide$rel_duration_B) / 2

# mean of both cohabiting duration reported
df_wide$hab_duration_A <- df_wide$pre_hab_duration_y_A + df_wide$pre_hab_duration_m_A / 12
df_wide$hab_duration_B <- df_wide$pre_hab_duration_y_B + df_wide$pre_hab_duration_m_B / 12 
df_wide$Cohabiting_Duration <- (df_wide$hab_duration_A + df_wide$hab_duration_B) / 2

#Gender of Each Partner
df_wide$Gender_A <- factor(df_wide$gender_A, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
df_wide$Gender_B <- factor(df_wide$gender_B, levels = c(1,2,3), labels = c('Male','Female', 'Other'))
#Same sex couples or not
df_wide$Same_Sex <- factor(df_wide$Gender_A == df_wide$Gender_B, levels = c(FALSE, TRUE), labels = c('Same-Sex Couple', 'Mixed-Sex Couple'))

#Handedness
df_wide$Handedness_A <- factor(df_wide$pre_handedness_A, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))
df_wide$Handedness_B <- factor(df_wide$pre_handedness_B, levels = c(0, 1, 2), c('Right','Left', 'Ambidextrous'))



# Income (mean of both partner's report)

merge_income <- function(income1, income2) {
  merged_income <- numeric(length(income1))
  
  # Loop through each pair of incomes
  for (i in seq_along(income1)) {
    # Handle NA
    if (is.na(income1[i])) {
      merged_income[i] <- income2[i]
    } 
    
    else if (is.na(income2[i])) {
      merged_income[i] <- income1[i]
    }
    
    # if both are informative, take mean and round
    else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
      merged_income[i] <- round((income1[i] + income2[i]) / 2)
    }
    
    # if one is informative and the other not, use the informative one. 
    else if (income1[i] %in% 1:6) {
      merged_income[i] <- income1[i]
    }
    else if (income2[i] %in% 1:6) {
      merged_income[i] <- income2[i]
    }
    
    # Now we only have cases, where both are either 70 or 99. We simply report "undisclosed". 
    else {
      merged_income[i] <- 99
    }
  }
  
  return(merged_income)
}


# Apply the function
numeric_income <- merge_income(df_wide$pre_income_1_A, df_wide$pre_income_1_B)

# Convert to factor
df_wide$Household_Income <- factor(numeric_income, levels = c(1,2,3,4,5,6,70,99), labels = c(
    "up to CHF 2'000.-", 
    "CHF 2'001.- to CHF 4'000.-",
    "CHF 4'001.- to CHF 6'000.-",
    "CHF 6'001.- to CHF 8'000.-",
    "CHF 8'001.- to CHF 10'000.-", 
    "above CHF 10'000.-", 
    "I don't know",
    "Undisclosed"
))





# Education for both separate
df_wide$Highest_Education_A <- factor(df_wide$pre_education_A, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)

df_wide$Highest_Education_B <- factor(df_wide$pre_education_B, levels = c(1,2,3,4,5,6,7), labels = c(
  "(still) no school diploma",
  "compulsory education (9 years)",
  "vocational training (apprenticeship)",
  "Matura (university entrance qualification)",
  "Bachelor's degree", 
  "Master's degree",
  "Doctorate degree"
  )
)



# Marital Status (married)
df_wide$Marital_Status_A <- df_wide$pre_mat_stat_A == 1
df_wide$Marital_status_B <- df_wide$pre_mat_stat_B == 1

# if at least one of them said they are married, we go for married
df_wide$Marital_Status <- factor((df_wide$Marital_Status_A + df_wide$Marital_status_B) > 0, levels = c(FALSE, TRUE), labels = c('Not Married', 'Married'))



# Children
df_wide$Children <- factor((df_wide$pre_child_option_A + df_wide$pre_child_option_B) > 0, levels = c(FALSE, TRUE), labels = c('Have Children', 'No Children'))
# number of chilren (ONLY FOR COUPLES THAT HAVE CHILDREN)
df_wide$nChildren_couples_with_children <- (df_wide$pre_child_nr_childs_A + df_wide$pre_child_nr_childs_B) / 2 
# number of children if every couple is considered. 
df_wide$nChildren_all_couples <- df_wide$nChildren_couples_with_children
df_wide$nChildren_all_couples[is.na(df_wide$nChildren_all_couples)] <- 0


# BMI (kg/m^2)
df_wide$pre_height_A <- df_wide$pre_height_A / 100 # to meters
df_wide$pre_height_B <- df_wide$pre_height_B / 100 # to meters

df_wide$BMI_A <- df_wide$pre_weight_A / (df_wide$pre_height_A^2)
df_wide$BMI_B <- df_wide$pre_weight_B / (df_wide$pre_height_B^2)



sample_measures <- c(
  "Gender_A", "Gender_B", "Same_Sex",
  "pre_age_A", "pre_age_B", "BMI_A", "BMI_B", "Handedness_A", "Handedness_B", # physical stats
  "Household_Income", "Highest_Education_A", "Highest_Education_B", # SES
  "Marital_Status", "Relationship_Duration", "Cohabiting_Duration",
  "Children", "nChildren_all_couples", "nChildren_couples_with_children" 
  )



sample_table <- report_measures(df_wide, sample_measures, ICC=F)


rownames(sample_table) <- NULL

openxlsx::write.xlsx(sample_table, 'Output/AsPreregistred/Sample.xlsx')

print_df(sample_table)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
Gender_A Other 0 0% NA NA NA NA
Gender_A Female 825 39% NA NA NA NA
Gender_A Male 1265 61% NA NA NA NA
Gender_B Other 0 0% NA NA NA NA
Gender_B Male 825 39% NA NA NA NA
Gender_B Female 1265 61% NA NA NA NA
Same_Sex Mixed-Sex Couple 110 5% NA NA NA NA
Same_Sex Same-Sex Couple 1980 95% NA NA NA NA
pre_age_A NA 2090 NA 0% 34.05 11.19 19.00-58.00
pre_age_B NA 2090 NA 0% 33.97 10.72 19.00-60.00
BMI_A NA 2090 NA 0% 24.58 3.77 16.37-33.95
BMI_B NA 2090 NA 0% 25.29 4.33 17.32-33.95
Handedness_A Ambidextrous 55 3% NA NA NA NA
Handedness_A Left 330 16% NA NA NA NA
Handedness_A Right 1705 82% NA NA NA NA
Handedness_B Ambidextrous 0 0% NA NA NA NA
Handedness_B Left 275 13% NA NA NA NA
Handedness_B Right 1815 87% NA NA NA NA
Household_Income I don’t know 0 0% NA NA NA NA
Household_Income up to CHF 2’000.- 110 5% NA NA NA NA
Household_Income CHF 2’001.- to CHF 4’000.- 165 8% NA NA NA NA
Household_Income CHF 6’001.- to CHF 8’000.- 165 8% NA NA NA NA
Household_Income CHF 4’001.- to CHF 6’000.- 275 13% NA NA NA NA
Household_Income Undisclosed 275 13% NA NA NA NA
Household_Income CHF 8’001.- to CHF 10’000.- 440 21% NA NA NA NA
Household_Income above CHF 10’000.- 660 32% NA NA NA NA
Highest_Education_A (still) no school diploma 0 0% NA NA NA NA
Highest_Education_A compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education_A Doctorate degree 0 0% NA NA NA NA
Highest_Education_A vocational training (apprenticeship) 330 16% NA NA NA NA
Highest_Education_A Master’s degree 495 24% NA NA NA NA
Highest_Education_A Matura (university entrance qualification) 605 29% NA NA NA NA
Highest_Education_A Bachelor’s degree 660 32% NA NA NA NA
Highest_Education_B (still) no school diploma 0 0% NA NA NA NA
Highest_Education_B compulsory education (9 years) 0 0% NA NA NA NA
Highest_Education_B Doctorate degree 0 0% NA NA NA NA
Highest_Education_B vocational training (apprenticeship) 110 5% NA NA NA NA
Highest_Education_B Bachelor’s degree 605 29% NA NA NA NA
Highest_Education_B Master’s degree 660 32% NA NA NA NA
Highest_Education_B Matura (university entrance qualification) 715 34% NA NA NA NA
Marital_Status Married 715 34% NA NA NA NA
Marital_Status Not Married 1375 66% NA NA NA NA
Relationship_Duration NA 2090 NA 0% 9.23 9.03 0.58-36.00
Cohabiting_Duration NA 2090 NA 0% 7.53 9.14 0.25-33.00
Children No Children 605 29% NA NA NA NA
Children Have Children 1485 71% NA NA NA NA
nChildren_all_couples NA 2090 NA 0% 0.47 0.85 0.00- 3.00
nChildren_couples_with_children NA 2090 NA 74% 1.80 0.60 1.00- 3.00

Descriptives of Measures

Main Measures

# Main Measures
main_constructs <- c("persuasion_A", "persuasion_B", 
                     "pressure_A", "pressure_B", 
                     "pushing_A", "pushing_B",
                     
                     "aff_A", "aff_B", 
                     "pa_sub_A", "pa_sub_B","pa_obj_A", "pa_obj_B", 
                     "anticipAff_A", "anticipAff_B", 
                     "reactance_A", "reactance_B"
                     )

main_descriptives <- report_measures(df_wide, main_constructs)

openxlsx::write.xlsx(main_descriptives, 'Output/AsPreregistred/DescriptivesMain.xlsx')

print_df(main_descriptives)
Variable n_Obs Missing Mean SD Range ICC
6 persuasion_A 2090 7% 0.43 1.09 0.00- 5.00 0.24
5 persuasion_B 2090 5% 0.41 1.08 0.00- 5.00 0.23
2 pressure_A 2090 7% 0.13 0.62 0.00- 5.00 0.63
1 pressure_B 2090 5% 0.12 0.62 0.00- 5.00 0.53
4 pushing_A 2090 7% 0.18 0.68 0.00- 5.00 0.09
3 pushing_B 2090 5% 0.15 0.64 0.00- 5.00 0.12
12 aff_A 2090 6% 4.89 1.13 1.00- 6.00 0.47
11 aff_B 2090 5% 4.77 1.14 1.00- 6.00 0.36
14 pa_sub_A 2090 6% 31.66 56.73 0.00-720.00 0.17
13 pa_sub_B 2090 5% 29.16 52.77 0.00-480.00 0.15
16 pa_obj_A 2090 12% 164.63 133.55 7.75-971.25 0.14
15 pa_obj_B 2090 11% 124.32 95.63 5.75-925.25 0.20
10 anticipAff_A 2090 6% 1.50 2.23 -5.00- 5.00 0.39
9 anticipAff_B 2090 5% 1.25 2.39 -5.00- 5.00 0.37
7 reactance_A 2090 80% 0.70 1.24 0.00- 4.00 0.52
8 reactance_B 2090 82% 0.90 1.39 0.00- 5.00 0.42

Including all Covariates

# With Covariates


# JITAIs as factor
df_wide$got_JITAI_A_factor <- factor(df_wide$got_JITAI_A, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))
df_wide$got_JITAI_B_factor <- factor(df_wide$got_JITAI_B, levels = c(0,1), labels = c("No JITAI", "Got a JITAI"))


# weekend as factor
df_wide$isWeekend_factor <- factor(df_wide$isWeekend_A, levels = c(0,1), labels = c("weekday", "weekend day"))


# Plan as a factor
df_wide$plan_A_factor <- factor(df_wide$plan_A, levels = c(0,1), labels = c("No Plan", "Yes Plan"))
df_wide$plan_B_factor <- factor(df_wide$plan_B, levels = c(0,1), labels = c("No Plan", "Yes Plan"))


# Study Group as a Factor
df_wide$studyGroup_factor <- factor(df_wide$studyGroup_A - 1, levels = c(0,1,2), labels = c("Allways Interventions", "First 3 weeks", "last 3 weeks"))

# Skilled support as factor
df_wide$skilled_support_factor <- factor(df_wide$skilled_support_A, levels = c(0,1), labels = c("Did not yet have SS", "Did already have SS"))



all_constructs <- c(
  main_constructs,
  "day",
  "weartime_A", "weartime_B",
  "barriers_A" , "barriers_B",
  "isWeekend_factor",
  "plan_A_factor", "plan_B_factor",
  "studyGroup_factor",
  "support_A", "support_B",
  "got_JITAI_A_factor", "got_JITAI_B_factor",
  "skilled_support_factor"
)


all_descriptives <- report_measures(df_wide, all_constructs, ICC = F)

openxlsx::write.xlsx(all_descriptives, 'Output/AsPreregistred/DescriptivesAll.xlsx')

print_df(all_descriptives)
Variable Level n_Obs percentage_Obs Missing Mean SD Range
29 persuasion_A NA 2090 NA 7% 0.43 1.09 0.00- 5.00
30 persuasion_B NA 2090 NA 5% 0.41 1.08 0.00- 5.00
31 pressure_A NA 2090 NA 7% 0.13 0.62 0.00- 5.00
32 pressure_B NA 2090 NA 5% 0.12 0.62 0.00- 5.00
33 pushing_A NA 2090 NA 7% 0.18 0.68 0.00- 5.00
34 pushing_B NA 2090 NA 5% 0.15 0.64 0.00- 5.00
18 aff_A NA 2090 NA 6% 4.89 1.13 1.00- 6.00
19 aff_B NA 2090 NA 5% 4.77 1.14 1.00- 6.00
27 pa_sub_A NA 2090 NA 6% 31.66 56.73 0.00- 720.00
28 pa_sub_B NA 2090 NA 5% 29.16 52.77 0.00- 480.00
25 pa_obj_A NA 2090 NA 12% 164.63 133.55 7.75- 971.25
26 pa_obj_B NA 2090 NA 11% 124.32 95.63 5.75- 925.25
20 anticipAff_A NA 2090 NA 6% 1.50 2.23 -5.00- 5.00
21 anticipAff_B NA 2090 NA 5% 1.25 2.39 -5.00- 5.00
35 reactance_A NA 2090 NA 80% 0.70 1.24 0.00- 4.00
36 reactance_B NA 2090 NA 82% 0.90 1.39 0.00- 5.00
24 day NA 2090 NA 0% 0.50 0.29 0.00- 1.00
39 weartime_A NA 2090 NA 2% 1036.85 374.95 0.00-1440.00
40 weartime_B NA 2090 NA 0% 1077.63 392.29 0.00-1440.00
22 barriers_A NA 2090 NA 54% 0.04 1.58 -3.29- 5.00
23 barriers_B NA 2090 NA 57% 0.10 1.70 -3.57- 5.00
6 isWeekend_factor weekend day 608 29% NA NA NA NA
14 isWeekend_factor weekday 1482 71% NA NA NA NA
2 plan_A_factor missing 136 7% NA NA NA NA
11 plan_A_factor Yes Plan 959 46% NA NA NA NA
12 plan_A_factor No Plan 995 48% NA NA NA NA
1 plan_B_factor missing 102 5% NA NA NA NA
10 plan_B_factor Yes Plan 901 43% NA NA NA NA
13 plan_B_factor No Plan 1087 52% NA NA NA NA
7 studyGroup_factor last 3 weeks 660 32% NA NA NA NA
8 studyGroup_factor Allways Interventions 715 34% NA NA NA NA
9 studyGroup_factor First 3 weeks 715 34% NA NA NA NA
37 support_A NA 2090 NA 7% 0.87 1.45 0.00- 5.00
38 support_B NA 2090 NA 5% 0.94 1.58 0.00- 5.00
3 got_JITAI_A_factor Got a JITAI 287 14% NA NA NA NA
17 got_JITAI_A_factor No JITAI 1803 86% NA NA NA NA
4 got_JITAI_B_factor Got a JITAI 298 14% NA NA NA NA
16 got_JITAI_B_factor No JITAI 1792 86% NA NA NA NA
5 skilled_support_factor Did not yet have SS 518 25% NA NA NA NA
15 skilled_support_factor Did already have SS 1572 75% NA NA NA NA

Correlations Main Constructs

cors <- wbCorr(df_wide[,c(main_constructs)], df_wide$coupleID, method = 'spearman')

main_cors <- summary(cors, 'wb')$merged_wb


openxlsx::write.xlsx(main_cors, 'Output/AsPreregistred/Correlations.xlsx')

print_df(main_cors, width = '7em')
persuasion_A persuasion_B pressure_A pressure_B pushing_A pushing_B aff_A aff_B pa_sub_A pa_sub_B pa_obj_A pa_obj_B anticipAff_A anticipAff_B reactance_A reactance_B
persuasion_A [0.24] 0.23*** 0.29*** 0.19*** 0.38*** 0.07** -0.02 0.03 0.14*** 0.11*** 0.01 0.00 0.02 0.03 -0.08 -0.01
persuasion_B 0.32 [0.23] 0.03 0.27*** 0.08*** 0.37*** -0.01 -0.01 0.20*** 0.20*** 0.07** 0.13*** 0.05* 0.04 0.01 0.02
pressure_A 0.46** 0.02 [0.63] 0.30*** 0.33*** -0.00 0.01 0.01 -0.01 -0.01 -0.03 -0.01 0.07** 0.06** 0.24*** 0.03
pressure_B 0.47** 0.34* 0.45** [0.53] 0.05* 0.23*** -0.02 -0.02 0.02 -0.04* -0.08** 0.06* 0.04 0.02 0.05 0.28***
pushing_A 0.55*** 0.23 0.48** 0.24 [0.09] 0.17*** -0.01 0.02 0.10*** 0.10*** 0.02 -0.02 0.02 0.02 0.09 -0.11
pushing_B 0.24 0.67*** 0.11 0.38* 0.35* [0.12] 0.03 0.02 0.17*** 0.13*** 0.07** 0.10*** 0.04 0.04 0.01 0.00
aff_A 0.35* 0.19 -0.03 -0.06 0.24 -0.03 [0.47] 0.19*** 0.19*** 0.11*** 0.13*** 0.12*** 0.30*** 0.10*** -0.02 0.00
aff_B 0.23 -0.04 0.01 -0.12 0.35* -0.01 0.41* [0.36] 0.10*** 0.20*** 0.10*** 0.18*** 0.08*** 0.36*** -0.02 -0.04
pa_sub_A 0.12 0.23 -0.27 -0.06 0.13 0.20 0.54*** 0.47** [0.17] 0.56*** 0.28*** 0.24*** 0.19*** 0.11*** -0.09 0.02
pa_sub_B 0.08 0.26 -0.16 -0.19 0.13 0.11 0.19 0.46** 0.70*** [0.15] 0.22*** 0.33*** 0.11*** 0.20*** 0.04 -0.05
pa_obj_A 0.03 0.22 0.02 -0.22 0.26 0.20 0.12 0.25 0.30 0.38* [0.14] 0.23*** 0.05* 0.03 0.01 0.03
pa_obj_B -0.06 0.18 0.04 -0.20 0.19 0.10 -0.00 0.40* 0.37* 0.54*** 0.76*** [0.20] 0.11*** 0.11*** 0.11* 0.13*
anticipAff_A 0.24 0.40* -0.35* -0.06 0.14 0.20 0.47** 0.38* 0.66*** 0.49** 0.40* 0.39* [0.39] 0.15*** -0.05 0.04
anticipAff_B 0.44** 0.12 0.07 0.25 0.39* 0.17 0.10 0.62*** 0.29 0.44** 0.11 0.28 0.41** [0.37] 0.00 -0.01
reactance_A 0.07 0.03 0.32* 0.07 0.14 0.03 -0.16 0.03 -0.13 -0.09 0.40* 0.37* 0.03 0.04 [0.52] 0.26***
reactance_B 0.17 0.30 -0.01 0.22 -0.16 0.29 -0.04 -0.15 -0.01 -0.03 0.07 0.06 0.11 0.02 0.01 [0.42]
plot(cors,'w')
## This may take a while...

#Modelling

Writing custion summary functions to facilitate reporting

summarize_brms <- function(model, short_version = FALSE, exponentiate = FALSE, model_rows_fixed = NULL, model_rows_random = NULL) {

  summ_og <- summary(model)
  summ <- summ_og$fixed
  rands <- summ_og$random[[1]][grep('sd\\(', rownames(summ_og$random[[1]])), ]
  
  ar1sterr <- summ_og$cor_pars
  sdresid <- summ_og$spec_pars
  
  rands <- rbind(rands, ar1sterr, sdresid)
  
  # Get rows that we need
  if (is.null(model_rows_fixed)) {
    model_rows_fixed <- rownames(summ)
  }
  if (is.null(model_rows_random)) {
    model_rows_random <- rownames(rands)
  }
  
  summ <- summ[model_rows_fixed, ]
  rands <- format(round(rands[model_rows_random, ], 2), nsmall = 2)
  
  Rhat <- format(round(summ$Rhat, 3), nsmall = 3)

  summ$Est.Error <- format(round(summ$Est.Error, 2), nsmall=2)
  summ$Bulk_ESS <- format(round(summ$Bulk_ESS, 2), nsmall=2)
  summ$Tail_ESS <- format(round(summ$Tail_ESS, 2), nsmall=2)
  
  #exponentiate if needed
  if (exponentiate) {
    summ$Estimate <- exp(summ$Estimate)
  }
    
  # Add stars based on Credible Interval
  significance <- ifelse(
    (summ$`l-95% CI` > 0 & summ$`u-95% CI` > 0) | 
    (summ$`l-95% CI` < 0 & summ$`u-95% CI` < 0), 
    '*', 
    '')
  
  # exponentiate CI
  if (exponentiate) {
    summ$`u-95% CI`<- exp(summ$`u-95% CI`)
    summ$`l-95% CI` <- exp(summ$`l-95% CI`)
  }
      

  
  # Round CI
  summ$`l-95% CI` <- format(round(summ$`l-95% CI`, 2), nsmall = 2)
  summ$`u-95% CI` <- format(round(summ$`u-95% CI`, 2), nsmall = 2)
  
  
  summ$Estimate <- ifelse( 
    is.na(summ$Estimate), 
    NA, 
    paste0(
      format(round(summ$Estimate,2), nsmall=2), 
      significance
    )
  )
  
  
  # Rename "Estimate" column correctly
  if (exponentiate) {
    if ('bernoulli' %in% model$family) {
      correct_name <- 'OR'
    } else if ('negbinomial' %in% model$family) {
      correct_name <- 'IRR'
    } else {
      correct_name <- 'exp(Est.)'
      warning(
        "Coefficients were exponentiated. Double check if this was intended."
        )
    }
    
  } else {
    correct_name <- 'b'
  }
  
  colnames(summ)[1] <- correct_name
  colnames(rands)[1] <- correct_name
  
  
  # Add rhat rounded to 3 digits back
  summ$Rhat <- Rhat
  
  # Add back all rownames
  rownames(summ) <- model_rows_fixed
  rownames(rands) <- model_rows_random
  
  # Add random effects to fixed effects df. 
  fullrep <- rbind(summ, rands)
  

  
  # We are finished, if we want the full table
  if (!short_version) {
      return(fullrep[ , c(1, 3, 4, 5, 6, 7)])
  }
  
  
  fullrep$`95% CI` <-  ifelse(
    is.na(fullrep[correct_name]) | fullrep$`u-95% CI` == "NA",
    NA,
    paste0("[", fullrep$`l-95% CI`, ", ", fullrep$`u-95% CI`, "]")
  )

  return(fullrep[ , c(1, 8)])
}


report_side_by_side <- function(..., model_rows_random, model_rows_fixed) {
  models <- list(...)
  model_names <- sapply(substitute(list(...))[-1], deparse)
  names(models) <- model_names

  side_by_side <- NULL
  for (i in seq_along(models)) {
    model <- models[[i]]
    model_name <- model_names[i]
    print(model_name)
    if ('bernoulli' %in% model$family | 'negbinomial' %in% model$family | grepl('log', model_name)) {
      exponentiate <- TRUE
    } else {
      exponentiate <- FALSE
    }
    model_summary <- summarize_brms(
      model, short_version = TRUE, 
      exponentiate = exponentiate, 
      model_rows_random = model_rows_random,
      model_rows_fixed = model_rows_fixed
      )
    
    colnames(model_summary) <- paste(colnames(model_summary), model_name)
    
    if (is.null(side_by_side)) {
      side_by_side <- model_summary
    } else {
      side_by_side <- cbind(side_by_side, model_summary)
    }
  }
  return(side_by_side)
}

MODELLING Persuasion

# For indistinguishable Dyads
model_rows_fixed_persuasion_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_persuasion_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_persuasion_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:persuasion_self_cw',
    'is_B:persuasion_self_cw', 
    'is_A:persuasion_partner_cw', 
    'is_B:persuasion_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:persuasion_self_cb',
    'is_B:persuasion_self_cb', 
    'is_A:persuasion_partner_cb', 
    'is_B:persuasion_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_persuasion_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:persuasion_self_cw)',
  'sd(is_B:persuasion_self_cw)', 
  'sd(is_A:persuasion_partner_cw)', 
  'sd(is_B:persuasion_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ persuasion

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros, transformations won’t help estimating the models. We employ the negative binomial family.

Indistinguishable

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + 
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "persuasion_pa_sub_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_sub_ind)
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12082.9 177.5
## p_loo        28.2   2.4
## looic     24165.8 355.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_sub_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 29.14* 20.18 41.57 1.001 903.07 1765.89
Within-Person Effects
persuasion_self_cw 1.01 0.92 1.11 1.000 9103.37 3111.55
persuasion_partner_cw 1.08 0.98 1.19 1.000 7783.86 2766.72
day 1.15 0.85 1.57 1.002 9017.24 3175.00
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.85 0.60 1.18 1.001 4385.12 2696.25
persuasion_partner_cb 1.11 0.79 1.59 1.000 3990.53 3056.64
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.72 0.52 0.98 1.01 1069.15 2159.42
sd(persuasion_self_cw) 0.25 0.12 0.40 1.00 2049.24 2247.12
sd(persuasion_partner_cw) 0.23 0.10 0.38 1.00 1924.16 2229.31
Additional Parameters
ar[1] 0.02 -0.94 0.94 1.00 6715.11 2833.16
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.00 8147.96 2612.93
sderr 0.05 0.00 0.13 1.00 2816.69 2191.39
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_sub_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_sub_apim)
## Warning: Found 4 observations with a pareto_k > 0.7 in model 'persuasion_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12085.1 177.4
## p_loo        40.5   2.7
## looic     24170.2 354.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3738  99.9%   206     
##    (0.7, 1]   (bad)         4   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_sub_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 31.64* 21.06 47.91 1.001 1568.57 2509.33
is_B 20.52* 14.28 29.88 1.000 2136.80 2886.58
Within-Person Effects
is_A:persuasion_self_cw 0.98 0.86 1.14 1.002 7023.50 3392.76
is_B:persuasion_self_cw 1.05 0.91 1.21 1.001 6912.83 3178.03
is_A:persuasion_partner_cw 0.99 0.87 1.14 1.001 8085.20 2979.63
is_B:persuasion_partner_cw 1.17* 1.01 1.36 1.000 6556.89 3694.90
is_A:day 1.19 0.75 1.88 1.000 5424.27 2903.84
is_B:day 1.27 0.82 1.98 1.001 6204.54 3147.04
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.93 0.57 1.49 1.001 3673.30 3212.62
is_B:persuasion_self_cb 0.64 0.37 1.09 1.001 3619.34 3300.48
is_A:persuasion_partner_cb 1.09 0.65 1.87 1.000 3627.47 2952.14
is_B:persuasion_partner_cb 1.16 0.71 1.91 1.001 3721.53 3193.61
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.72 0.51 0.97 1.00 1950.14 2855.84
sd(is_B) 0.64 0.43 0.89 1.00 1661.64 2675.85
sd(is_A:persuasion_self_cw) 0.22 0.04 0.41 1.00 1392.37 1070.47
sd(is_B:persuasion_self_cw) 0.23 0.04 0.44 1.00 1545.44 1293.99
sd(is_A:persuasion_partner_cw) 0.18 0.01 0.40 1.00 1508.09 1470.95
sd(is_B:persuasion_partner_cw) 0.25 0.06 0.46 1.00 1662.67 1452.00
Additional Parameters
ar[1] 0.01 -0.94 0.96 1.00 5779.57 2085.73
nu NA NA NA NA NA NA
shape 0.14 0.13 0.14 1.00 6750.82 2677.34
sderr 0.04 0.00 0.13 1.00 2791.50 2417.53
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                        elpd_diff se_diff
## persuasion_pa_sub_ind   0.0       0.0   
## persuasion_pa_sub_apim -2.2       3.7
bayes_R2(persuasion_pa_sub_apim)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.2106389 0.0702662 0.1169052 0.4020931
bayes_R2(persuasion_pa_sub_ind)
##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1703761 0.05028759 0.09973377 0.2920029

Device Based MVPA ~ persuasion

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the model did not converge. Poisson also did not work. As we have no zeros in this distribution, we log transform.

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_obj_log_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_pa_obj_log_ind)
print(loo_ind)
## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2803.1  54.8
## p_loo        70.8   3.0
## looic      5606.2 109.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.5]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(persuasion_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_persuasion_ind, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 120.38* 108.19 134.27 1.011 458.70 1046.96
Within-Person Effects
persuasion_self_cw 1.01 0.99 1.03 1.001 6192.13 3094.03
persuasion_partner_cw 1.01 0.99 1.02 1.002 5939.04 3243.98
day 0.99 0.91 1.08 1.001 4877.50 3218.69
weartime_self_cw 1.00 1.00 1.00 1.002 3797.90 2901.27
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.87* 0.79 0.97 1.001 2456.65 2873.55
persuasion_partner_cb 1.05 0.95 1.17 1.001 2390.52 2787.94
weartime_self_cb 1.00 1.00 1.00 1.001 4350.93 3515.45
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 0.23 0.38 1.01 758.27 1077.55
sd(persuasion_self_cw) 0.06 0.03 0.09 1.00 2263.37 2401.22
sd(persuasion_partner_cw) 0.05 0.03 0.08 1.00 2075.45 1844.66
Additional Parameters
ar[1] 0.30 0.26 0.33 1.00 5720.76 2947.82
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.56 0.55 0.57 1.00 5237.61 3296.42

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_pa_obj_log_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_pa_obj_log_apim)
print(loo_apim)
## 
## Computed from 4000 by 3305 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2721.8  56.3
## p_loo       129.9   5.7
## looic      5443.6 112.5
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.2]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  persuasion_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(persuasion_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_persuasion_apim,
## : Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 129.47* 108.37 153.25 1.005 675.45 1343.39
is_B 118.68* 101.47 140.04 1.003 717.81 1147.66
Within-Person Effects
is_A:persuasion_self_cw 0.97 0.95 1.00 1.001 6903.31 2865.10
is_B:persuasion_self_cw 1.05* 1.02 1.07 1.000 8393.27 3023.29
is_A:persuasion_partner_cw 1.00 0.97 1.03 1.000 10350.62 2548.48
is_B:persuasion_partner_cw 1.01 0.98 1.04 1.001 8481.91 2454.68
is_A:day 0.89* 0.79 0.99 1.000 5695.72 2639.66
is_B:day 1.06 0.95 1.19 1.004 6999.61 3316.51
is_A:weartime_self_cw 1.00 1.00 1.00 1.000 4098.01 2327.61
is_B:weartime_self_cw 1.00 1.00 1.00 1.001 3896.57 3064.60
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.89 0.77 1.04 1.000 2996.19 3204.21
is_B:persuasion_self_cb 0.85 0.71 1.00 1.000 2505.04 2772.26
is_A:persuasion_partner_cb 1.08 0.93 1.26 1.001 3333.22 2973.30
is_B:persuasion_partner_cb 1.15 0.99 1.34 1.001 2837.75 2698.05
is_A:weartime_self_cb 1.00 1.00 1.00 1.000 4489.33 3458.59
is_B:weartime_self_cb 1.00* 1.00 1.00 1.000 4627.69 3711.40
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.31 0.23 0.42 1.00 918.75 1454.42
sd(is_B) 0.39 0.29 0.52 1.00 736.28 1291.88
sd(is_A:persuasion_self_cw) 0.05 0.01 0.10 1.00 1396.31 1172.67
sd(is_B:persuasion_self_cw) 0.08 0.05 0.13 1.00 2783.30 3043.30
sd(is_A:persuasion_partner_cw) 0.07 0.03 0.13 1.00 1588.95 1861.31
sd(is_B:persuasion_partner_cw) 0.05 0.00 0.10 1.00 1064.17 842.20
Additional Parameters
ar[1] 0.22 0.18 0.25 1.00 6900.91 2930.32
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.54 0.53 0.55 1.00 7892.53 2411.23

Compare Models

loo_compare(loo_ind, loo_apim)
##                            elpd_diff se_diff
## persuasion_pa_obj_log_apim   0.0       0.0  
## persuasion_pa_obj_log_ind  -81.3      14.3
bayes_R2(persuasion_pa_obj_log_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1890189 0.006237183 0.1766613 0.2013082
bayes_R2(persuasion_pa_obj_log_ind)
##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.1678555 0.006690044 0.1543916 0.180971

Affect ~ persuasion

range(df_double$aff, na.rm = T) 
## [1] 1 6
hist(df_double$aff, breaks = 15)

Indistinguishable

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_aff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_aff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4821.9  63.5
## p_loo        63.0   2.6
## looic      9643.9 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_aff_ind, ask = FALSE)

summarize_brms(
  persuasion_aff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.90* 4.70 5.11 1.002 760.07 1412.86
Within-Person Effects
persuasion_self_cw -0.01 -0.04 0.02 1.002 6135.29 2936.42
persuasion_partner_cw 0.03 0.00 0.06 1.001 7257.53 2884.58
day -0.10 -0.25 0.04 1.000 6739.82 3092.48
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.02 -0.16 0.21 1.002 2929.97 2540.56
persuasion_partner_cb 0.17 -0.02 0.35 1.003 2577.32 2745.95
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.48 0.77 1.00 1006.48 1542.43
sd(persuasion_self_cw) 0.03 0.00 0.07 1.00 1663.48 1453.32
sd(persuasion_partner_cw) 0.07 0.01 0.12 1.00 929.44 1025.45
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 6242.15 2817.43
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 5910.04 2965.81

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_aff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_aff_apim)
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4694.8  65.2
## p_loo       118.3   5.0
## looic      9389.6 130.5
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_aff_apim, ask = FALSE)

summarize_brms(
  persuasion_aff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 4.88* 4.61 5.14 1.001 713.17 1297.44
is_B 4.89* 4.61 5.16 1.001 771.84 1579.80
Within-Person Effects
is_A:persuasion_self_cw -0.02 -0.06 0.02 1.002 8436.28 2944.95
is_B:persuasion_self_cw 0.00 -0.04 0.04 1.000 9076.31 3247.40
is_A:persuasion_partner_cw 0.03 -0.01 0.07 1.003 7295.23 2741.16
is_B:persuasion_partner_cw 0.03 -0.01 0.07 1.002 7076.20 3073.96
is_A:day -0.03 -0.21 0.16 1.004 6493.63 2849.65
is_B:day -0.16 -0.34 0.02 1.003 8342.01 3166.59
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.05 -0.19 0.28 1.000 3092.29 3382.04
is_B:persuasion_self_cb -0.12 -0.41 0.17 1.001 2558.03 2660.06
is_A:persuasion_partner_cb 0.07 -0.18 0.32 1.001 3230.63 3103.19
is_B:persuasion_partner_cb 0.01 -0.26 0.29 1.002 2399.79 2699.47
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.80 0.63 1.03 1.00 1146.74 1999.75
sd(is_B) 0.68 0.53 0.87 1.00 955.61 1858.38
sd(is_A:persuasion_self_cw) 0.05 0.00 0.11 1.01 1036.51 1603.46
sd(is_B:persuasion_self_cw) 0.08 0.01 0.15 1.00 1336.77 1537.21
sd(is_A:persuasion_partner_cw) 0.05 0.00 0.12 1.00 1388.60 1697.67
sd(is_B:persuasion_partner_cw) 0.08 0.01 0.15 1.00 1193.25 1533.01
Additional Parameters
ar[1] 0.34 0.30 0.37 1.00 5874.13 2813.52
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.83 0.82 0.85 1.00 7624.11 3099.01

Compare Models

loo_compare(loo_ind, loo_apim)
##                     elpd_diff se_diff
## persuasion_aff_apim    0.0       0.0 
## persuasion_aff_ind  -127.1      18.9
bayes_R2(persuasion_aff_apim)
##     Estimate  Est.Error    Q2.5     Q97.5
## R2 0.2341884 0.00451868 0.22503 0.2427378
bayes_R2(persuasion_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2247693 0.004741361 0.2155174 0.2337927

Anticipatory Affect ~ persuasion

range(df_double$anticipAff, na.rm = T) 
## [1] -5  5
hist(df_double$anticipAff, breaks = 10)

Indistinguishable

formula <- bf(
  anticipAff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_anticipAff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7536.8  56.6
## p_loo        64.3   2.6
## looic     15073.6 113.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_anticipAff_ind, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.47* 1.01 1.93 1.012 374.59 740.14
Within-Person Effects
persuasion_self_cw -0.03 -0.08 0.03 1.000 4294.03 3090.91
persuasion_partner_cw 0.08* 0.02 0.14 1.000 4654.38 3167.04
day -0.01 -0.32 0.30 1.001 4311.50 2389.89
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb -0.21 -0.57 0.15 1.002 1726.57 2121.35
persuasion_partner_cb 0.43* 0.05 0.80 1.002 1930.37 2210.76
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.27 0.99 1.60 1.01 583.64 1095.04
sd(persuasion_self_cw) 0.16 0.07 0.27 1.00 1388.45 1511.56
sd(persuasion_partner_cw) 0.04 0.00 0.11 1.00 1559.41 1549.23
Additional Parameters
ar[1] 0.42 0.39 0.45 1.00 3932.52 2916.81
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.80 1.76 1.84 1.00 5348.52 2742.81

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 5000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_anticipAff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_anticipAff_apim)
print(loo_apim)
## 
## Computed from 4000 by 3742 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7439.2  56.7
## p_loo       114.1   4.2
## looic     14878.4 113.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_anticipAff_apim, ask = FALSE)

summarize_brms(
  persuasion_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.23* 0.60 1.86 1.004 607.10 1139.00
is_B 1.69* 1.15 2.27 1.002 817.88 1271.71
Within-Person Effects
is_A:persuasion_self_cw -0.01 -0.09 0.07 1.003 8297.26 2853.71
is_B:persuasion_self_cw -0.03 -0.11 0.05 1.003 8454.33 2857.25
is_A:persuasion_partner_cw 0.09* 0.01 0.17 1.001 8564.71 2885.68
is_B:persuasion_partner_cw 0.07 -0.01 0.15 1.000 8440.49 3254.58
is_A:day 0.06 -0.33 0.44 1.002 7256.56 2746.37
is_B:day -0.08 -0.45 0.30 1.002 6818.28 3537.89
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb -0.22 -0.70 0.27 1.001 3493.49 3307.53
is_B:persuasion_self_cb -0.59 -1.17 0.00 1.000 2676.13 2852.30
is_A:persuasion_partner_cb 0.43 -0.10 0.94 1.001 3815.72 3570.18
is_B:persuasion_partner_cb 0.50 -0.03 1.04 1.000 2770.49 2930.67
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 1.46 1.15 1.85 1.00 866.25 1926.05
sd(is_B) 1.49 1.16 1.91 1.00 1079.48 1532.73
sd(is_A:persuasion_self_cw) 0.13 0.01 0.29 1.00 1122.08 2193.39
sd(is_B:persuasion_self_cw) 0.20 0.07 0.34 1.00 1629.01 1530.04
sd(is_A:persuasion_partner_cw) 0.06 0.00 0.18 1.00 1975.80 2001.79
sd(is_B:persuasion_partner_cw) 0.06 0.00 0.17 1.00 2245.16 2141.55
Additional Parameters
ar[1] 0.33 0.30 0.37 1.00 7088.02 2876.23
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.74 1.70 1.78 1.00 8090.93 3063.62

Compare Models

loo_compare(loo_ind, loo_apim)
##                            elpd_diff se_diff
## persuasion_anticipAff_apim   0.0       0.0  
## persuasion_anticipAff_ind  -97.6      15.8
bayes_R2(persuasion_anticipAff_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2193938 0.004568784 0.2103348 0.2284799
bayes_R2(persuasion_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2105418 0.004928154 0.2010232 0.2202312

reactance ~ persuasion

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

Indistinguishable

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    persuasion_self_cb + persuasion_partner_cb +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_reactance_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(persuasion_reactance_ind)
print(loo_ind)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1101.1 33.5
## p_loo        54.8  5.3
## looic      2202.1 66.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(persuasion_reactance_ind, ask = FALSE)

summarize_brms(
  persuasion_reactance_ind, 
  model_rows_fixed = model_rows_fixed_persuasion_ind,
  model_rows_random = model_rows_random_persuasion_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.71* 0.43 0.98 1.002 1841.18 2594.40
Within-Person Effects
persuasion_self_cw -0.01 -0.06 0.05 1.001 7277.61 2455.13
persuasion_partner_cw -0.04 -0.10 0.02 1.002 9742.45 2931.45
day -0.15 -0.45 0.13 1.001 7404.57 3125.69
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.02 -0.21 0.25 1.000 4189.26 2960.39
persuasion_partner_cb 0.07 -0.18 0.31 1.000 4112.40 2982.04
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.46 0.79 1.00 1971.40 2825.33
sd(persuasion_self_cw) 0.11 0.02 0.22 1.00 671.59 973.55
sd(persuasion_partner_cw) 0.05 0.00 0.12 1.00 1730.57 1918.83
Additional Parameters
ar[1] 0.07 -0.01 0.14 1.00 6223.85 2979.48
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.00 0.95 1.06 1.00 4549.32 2608.52

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw +
    
    is_A:persuasion_self_cb + is_B:persuasion_self_cb + 
    is_A:persuasion_partner_cb + is_B:persuasion_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:persuasion_self_cw + is_B:persuasion_self_cw + 
    is_A:persuasion_partner_cw + is_B:persuasion_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

persuasion_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_persuasion_reactance_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(persuasion_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(persuasion_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(persuasion_reactance_apim)
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'persuasion_reactance_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1082.0 33.1
## p_loo        90.5  8.1
## looic      2164.1 66.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     753   99.6%   83      
##    (0.7, 1]   (bad)        3    0.4%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(persuasion_reactance_apim, ask = FALSE)

summarize_brms(
  persuasion_reactance_apim, 
  model_rows_fixed = model_rows_fixed_persuasion_apim,
  model_rows_random = model_rows_random_persuasion_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.40* 0.04 0.75 1.001 1679.69 2415.79
is_B 0.95* 0.58 1.33 1.001 2051.50 2436.77
Within-Person Effects
is_A:persuasion_self_cw 0.00 -0.08 0.07 1.000 4755.55 2853.10
is_B:persuasion_self_cw -0.01 -0.08 0.06 1.002 6490.63 2909.35
is_A:persuasion_partner_cw -0.02 -0.09 0.06 1.000 6363.01 3371.58
is_B:persuasion_partner_cw -0.08 -0.16 0.01 1.001 5639.58 2773.14
is_A:day 0.18 -0.21 0.57 1.000 4082.36 3290.51
is_B:day -0.37 -0.78 0.04 1.000 4530.83 3189.96
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.42* 0.09 0.76 1.001 2517.84 2853.79
is_B:persuasion_self_cb -0.26 -0.67 0.12 1.002 2994.88 2948.22
is_A:persuasion_partner_cb -0.53* -0.94 -0.14 1.001 2581.08 2864.87
is_B:persuasion_partner_cb 0.35* 0.03 0.68 1.001 2781.56 2553.57
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.70 0.51 0.96 1.00 1408.17 2080.97
sd(is_B) 0.66 0.49 0.91 1.00 1783.02 2768.51
sd(is_A:persuasion_self_cw) 0.21 0.07 0.35 1.00 680.11 868.96
sd(is_B:persuasion_self_cw) 0.08 0.00 0.19 1.01 1029.42 1675.40
sd(is_A:persuasion_partner_cw) 0.07 0.00 0.18 1.00 1891.93 2446.46
sd(is_B:persuasion_partner_cw) 0.10 0.01 0.22 1.00 1083.22 945.22
Additional Parameters
ar[1] -0.03 -0.12 0.05 1.00 3084.44 2569.49
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.95 0.89 1.00 1.00 3168.19 2657.30

Compare Models

loo_compare(loo_ind, loo_apim)
##                           elpd_diff se_diff
## persuasion_reactance_apim   0.0       0.0  
## persuasion_reactance_ind  -19.0       7.2
bayes_R2(persuasion_reactance_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2386339 0.01018588 0.2179162 0.2579992
bayes_R2(persuasion_reactance_ind)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.2194376 0.0112548 0.1976425 0.2413073

Report models for Persuasion side by side

Indistinguishable Models

persuasion_all_ind_models <- report_side_by_side(
  persuasion_pa_sub_ind,
  persuasion_pa_obj_log_ind,
  persuasion_aff_ind,
  persuasion_anticipAff_ind,
  persuasion_reactance_ind,
  
  model_rows_random = model_rows_random_persuasion_ind,
  model_rows_fixed = model_rows_fixed_persuasion_ind
) 
## [1] "persuasion_pa_sub_ind"
## [1] "persuasion_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "persuasion_aff_ind"
## [1] "persuasion_anticipAff_ind"
## [1] "persuasion_reactance_ind"
# pretty printing
persuasion_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR persuasion_pa_sub_ind 95% CI persuasion_pa_sub_ind exp(Est.) persuasion_pa_obj_log_ind 95% CI persuasion_pa_obj_log_ind b persuasion_aff_ind 95% CI persuasion_aff_ind b persuasion_anticipAff_ind 95% CI persuasion_anticipAff_ind b persuasion_reactance_ind 95% CI persuasion_reactance_ind
Intercept 29.14* [20.18, 41.57] 120.38* [108.19, 134.27] 4.90* [ 4.70, 5.11] 1.47* [ 1.01, 1.93] 0.71* [ 0.43, 0.98]
Within-Person Effects
persuasion_self_cw 1.01 [ 0.92, 1.11] 1.01 [ 0.99, 1.03] -0.01 [-0.04, 0.02] -0.03 [-0.08, 0.03] -0.01 [-0.06, 0.05]
persuasion_partner_cw 1.08 [ 0.98, 1.19] 1.01 [ 0.99, 1.02] 0.03 [ 0.00, 0.06] 0.08* [ 0.02, 0.14] -0.04 [-0.10, 0.02]
day 1.15 [ 0.85, 1.57] 0.99 [ 0.91, 1.08] -0.10 [-0.25, 0.04] -0.01 [-0.32, 0.30] -0.15 [-0.45, 0.13]
weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
persuasion_self_cb 0.85 [ 0.60, 1.18] 0.87* [ 0.79, 0.97] 0.02 [-0.16, 0.21] -0.21 [-0.57, 0.15] 0.02 [-0.21, 0.25]
persuasion_partner_cb 1.11 [ 0.79, 1.59] 1.05 [ 0.95, 1.17] 0.17 [-0.02, 0.35] 0.43* [ 0.05, 0.80] 0.07 [-0.18, 0.31]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.72 [ 0.52, 0.98] 0.30 [0.23, 0.38] 0.60 [0.48, 0.77] 1.27 [0.99, 1.60] 0.61 [ 0.46, 0.79]
sd(persuasion_self_cw) 0.25 [ 0.12, 0.40] 0.06 [0.03, 0.09] 0.03 [0.00, 0.07] 0.16 [0.07, 0.27] 0.11 [ 0.02, 0.22]
sd(persuasion_partner_cw) 0.23 [ 0.10, 0.38] 0.05 [0.03, 0.08] 0.07 [0.01, 0.12] 0.04 [0.00, 0.11] 0.05 [ 0.00, 0.12]
Additional Parameters
ar[1] 0.02 [-0.94, 0.94] 0.30 [0.26, 0.33] 0.45 [0.42, 0.48] 0.42 [0.39, 0.45] 0.07 [-0.01, 0.14]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.14 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.56 [0.55, 0.57] 0.87 [0.85, 0.89] 1.80 [1.76, 1.84] 1.00 [ 0.95, 1.06]
# prepare export to excel
all_models <- list()
all_models$Persuasion_ind <- prepare_for_export(persuasion_all_ind_models)

APIMs

persuasion_all_apim_models <- report_side_by_side(
  persuasion_pa_sub_apim,
  persuasion_pa_obj_log_apim,
  persuasion_aff_apim,
  persuasion_anticipAff_apim,
  persuasion_reactance_apim,
  
  model_rows_random = model_rows_random_persuasion_apim,
  model_rows_fixed = model_rows_fixed_persuasion_apim
) 
## [1] "persuasion_pa_sub_apim"
## [1] "persuasion_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "persuasion_aff_apim"
## [1] "persuasion_anticipAff_apim"
## [1] "persuasion_reactance_apim"
# pretty printing
persuasion_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR persuasion_pa_sub_apim 95% CI persuasion_pa_sub_apim exp(Est.) persuasion_pa_obj_log_apim 95% CI persuasion_pa_obj_log_apim b persuasion_aff_apim 95% CI persuasion_aff_apim b persuasion_anticipAff_apim 95% CI persuasion_anticipAff_apim b persuasion_reactance_apim 95% CI persuasion_reactance_apim
is_A 31.64* [21.06, 47.91] 129.47* [108.37, 153.25] 4.88* [ 4.61, 5.14] 1.23* [ 0.60, 1.86] 0.40* [ 0.04, 0.75]
is_B 20.52* [14.28, 29.88] 118.68* [101.47, 140.04] 4.89* [ 4.61, 5.16] 1.69* [ 1.15, 2.27] 0.95* [ 0.58, 1.33]
Within-Person Effects
is_A:persuasion_self_cw 0.98 [ 0.86, 1.14] 0.97 [ 0.95, 1.00] -0.02 [-0.06, 0.02] -0.01 [-0.09, 0.07] 0.00 [-0.08, 0.07]
is_B:persuasion_self_cw 1.05 [ 0.91, 1.21] 1.05* [ 1.02, 1.07] 0.00 [-0.04, 0.04] -0.03 [-0.11, 0.05] -0.01 [-0.08, 0.06]
is_A:persuasion_partner_cw 0.99 [ 0.87, 1.14] 1.00 [ 0.97, 1.03] 0.03 [-0.01, 0.07] 0.09* [ 0.01, 0.17] -0.02 [-0.09, 0.06]
is_B:persuasion_partner_cw 1.17* [ 1.01, 1.36] 1.01 [ 0.98, 1.04] 0.03 [-0.01, 0.07] 0.07 [-0.01, 0.15] -0.08 [-0.16, 0.01]
is_A:day 1.19 [ 0.75, 1.88] 0.89* [ 0.79, 0.99] -0.03 [-0.21, 0.16] 0.06 [-0.33, 0.44] 0.18 [-0.21, 0.57]
is_B:day 1.27 [ 0.82, 1.98] 1.06 [ 0.95, 1.19] -0.16 [-0.34, 0.02] -0.08 [-0.45, 0.30] -0.37 [-0.78, 0.04]
is_A:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
is_A:persuasion_self_cb 0.93 [ 0.57, 1.49] 0.89 [ 0.77, 1.04] 0.05 [-0.19, 0.28] -0.22 [-0.70, 0.27] 0.42* [ 0.09, 0.76]
is_B:persuasion_self_cb 0.64 [ 0.37, 1.09] 0.85 [ 0.71, 1.00] -0.12 [-0.41, 0.17] -0.59 [-1.17, 0.00] -0.26 [-0.67, 0.12]
is_A:persuasion_partner_cb 1.09 [ 0.65, 1.87] 1.08 [ 0.93, 1.26] 0.07 [-0.18, 0.32] 0.43 [-0.10, 0.94] -0.53* [-0.94, -0.14]
is_B:persuasion_partner_cb 1.16 [ 0.71, 1.91] 1.15 [ 0.99, 1.34] 0.01 [-0.26, 0.29] 0.50 [-0.03, 1.04] 0.35* [ 0.03, 0.68]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(is_A) 0.72 [ 0.51, 0.97] 0.31 [0.23, 0.42] 0.80 [0.63, 1.03] 1.46 [1.15, 1.85] 0.70 [ 0.51, 0.96]
sd(is_B) 0.64 [ 0.43, 0.89] 0.39 [0.29, 0.52] 0.68 [0.53, 0.87] 1.49 [1.16, 1.91] 0.66 [ 0.49, 0.91]
sd(is_A:persuasion_self_cw) 0.22 [ 0.04, 0.41] 0.05 [0.01, 0.10] 0.05 [0.00, 0.11] 0.13 [0.01, 0.29] 0.21 [ 0.07, 0.35]
sd(is_B:persuasion_self_cw) 0.23 [ 0.04, 0.44] 0.08 [0.05, 0.13] 0.08 [0.01, 0.15] 0.20 [0.07, 0.34] 0.08 [ 0.00, 0.19]
sd(is_A:persuasion_partner_cw) 0.18 [ 0.01, 0.40] 0.07 [0.03, 0.13] 0.05 [0.00, 0.12] 0.06 [0.00, 0.18] 0.07 [ 0.00, 0.18]
sd(is_B:persuasion_partner_cw) 0.25 [ 0.06, 0.46] 0.05 [0.00, 0.10] 0.08 [0.01, 0.15] 0.06 [0.00, 0.17] 0.10 [ 0.01, 0.22]
Additional Parameters
ar[1] 0.01 [-0.94, 0.96] 0.22 [0.18, 0.25] 0.34 [0.30, 0.37] 0.33 [0.30, 0.37] -0.03 [-0.12, 0.05]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.14 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.04 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.54 [0.53, 0.55] 0.83 [0.82, 0.85] 1.74 [1.70, 1.78] 0.95 [ 0.89, 1.00]
# prepare export to excel
all_models$Persuasion_apim <- prepare_for_export(persuasion_all_apim_models)

Modelling Pressure

# For indistinguishable Dyads
model_rows_fixed_pressure_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pressure_self_cb',
    'pressure_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pressure_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pressure_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pressure_self_cw',
    'is_B:pressure_self_cw', 
    'is_A:pressure_partner_cw', 
    'is_B:pressure_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pressure_self_cb',
    'is_B:pressure_self_cb', 
    'is_A:pressure_partner_cb', 
    'is_B:pressure_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pressure_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pressure_self_cw)',
  'sd(is_B:pressure_self_cw)', 
  'sd(is_A:pressure_partner_cw)', 
  'sd(is_B:pressure_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + 
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_sub_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_sub_ind)
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12087.0 177.4
## p_loo        20.0   1.5
## looic     24173.9 354.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 2.5]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_sub_ind, ask = FALSE)

summarize_brms(
  pressure_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 24.44* 18.79 32.14 1.001 2684.70 3008.73
Within-Person Effects
pressure_self_cw 0.98 0.79 1.23 1.001 9687.42 2765.59
pressure_partner_cw 0.86 0.69 1.10 1.003 9495.13 2763.51
day 1.18 0.86 1.62 1.002 9666.04 2220.42
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 1.09 0.66 1.81 1.000 5155.30 2888.05
pressure_partner_cb 0.89 0.55 1.47 1.000 5460.05 3143.76
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 0.47 0.84 1.00 1523.48 2255.08
sd(pressure_self_cw) 0.13 0.01 0.39 1.00 2670.03 2192.36
sd(pressure_partner_cw) 0.16 0.00 0.47 1.00 2200.86 1745.79
Additional Parameters
ar[1] 0.01 -0.93 0.94 1.00 7307.17 2816.12
nu NA NA NA NA NA NA
shape 0.13 0.12 0.14 1.00 9890.25 2099.07
sderr 0.04 0.00 0.13 1.00 2367.18 1497.39
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_sub_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_sub_apim)
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'pressure_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -12086.2 177.4
## p_loo        29.0   2.0
## looic     24172.4 354.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3731  99.9%   476     
##    (0.7, 1]   (bad)         5   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_sub_apim, ask = FALSE)

summarize_brms(
  pressure_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 29.61* 20.61 41.74 1.003 1776.95 2047.69
is_B 19.98* 14.03 28.42 1.002 1672.04 2176.58
Within-Person Effects
is_A:pressure_self_cw 0.78 0.59 1.07 1.001 6725.39 3351.31
is_B:pressure_self_cw 1.26 0.91 1.88 1.001 3866.85 1938.48
is_A:pressure_partner_cw 1.02 0.76 1.45 1.001 5490.25 2781.88
is_B:pressure_partner_cw 0.63* 0.44 0.95 1.000 4440.11 2881.48
is_A:day 1.07 0.69 1.68 1.000 4162.17 2860.47
is_B:day 1.20 0.78 1.81 1.003 4004.45 2855.89
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 1.08 0.49 2.36 1.000 1977.84 2075.41
is_B:pressure_self_cb 1.16 0.45 3.08 1.000 2129.31 2374.48
is_A:pressure_partner_cb 0.79 0.29 2.13 1.001 2090.44 2171.75
is_B:pressure_partner_cb 0.84 0.41 1.76 1.001 2526.13 2873.82
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.66 0.47 0.90 1.00 1972.83 2833.96
sd(is_B) 0.57 0.37 0.80 1.00 1870.38 2642.32
sd(is_A:pressure_self_cw) 0.21 0.01 0.70 1.00 2502.52 2167.92
sd(is_B:pressure_self_cw) 0.21 0.01 0.72 1.00 2768.34 2212.50
sd(is_A:pressure_partner_cw) 0.25 0.01 0.85 1.00 2796.74 2566.93
sd(is_B:pressure_partner_cw) 0.24 0.01 0.82 1.00 1578.33 1954.60
Additional Parameters
ar[1] 0.02 -0.94 0.94 1.00 4840.27 2464.21
nu NA NA NA NA NA NA
shape 0.13 0.13 0.14 1.00 4994.16 2779.02
sderr 0.04 0.00 0.12 1.00 2265.82 1795.20
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                      elpd_diff se_diff
## pressure_pa_sub_apim  0.0       0.0   
## pressure_pa_sub_ind  -0.7       2.6
bayes_R2(pressure_pa_sub_apim)
##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1782304 0.09020382 0.09058281 0.4855923
bayes_R2(pressure_pa_sub_ind)
##     Estimate  Est.Error       Q2.5     Q97.5
## R2 0.1167618 0.03779368 0.06777029 0.1964773

Device Based MVPA ~ pressure

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_obj_log_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_pa_obj_log_ind)
print(loo_ind)
## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2821.5  54.8
## p_loo        52.9   2.4
## looic      5643.0 109.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(pressure_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_pressure_ind, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 118.46* 106.30 131.66 1.009 535.38 1156.50
Within-Person Effects
pressure_self_cw 1.00 0.95 1.05 1.000 6377.09 3006.92
pressure_partner_cw 0.96 0.91 1.02 1.001 6449.17 2850.72
day 0.97 0.89 1.06 1.002 5721.89 2678.48
weartime_self_cw 1.00 1.00 1.00 1.001 3885.21 2144.43
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.78* 0.67 0.93 1.002 2469.56 2730.16
pressure_partner_cb 1.18* 1.02 1.37 1.002 2278.16 2571.02
weartime_self_cb 1.00 1.00 1.00 1.001 3872.79 3061.44
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.29 0.22 0.37 1.01 814.16 1512.49
sd(pressure_self_cw) 0.05 0.00 0.14 1.00 1496.27 1841.42
sd(pressure_partner_cw) 0.03 0.00 0.10 1.00 2216.36 2241.51
Additional Parameters
ar[1] 0.28 0.25 0.32 1.00 5765.91 2848.84
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.56 0.55 0.58 1.00 6556.32 3074.86

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_pa_obj_log_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_pa_obj_log_apim)
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'pressure_pa_obj_log_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3299 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2748.8  56.1
## p_loo       104.1   5.3
## looic      5497.7 112.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3294  99.8%   481     
##    (0.7, 1]   (bad)         5   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pressure_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(pressure_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_pressure_apim, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 127.96* 106.77 150.24 1.017 531.49 857.08
is_B 123.39* 104.27 145.22 1.012 424.74 1029.31
Within-Person Effects
is_A:pressure_self_cw 1.01 0.94 1.09 1.000 7896.70 2787.06
is_B:pressure_self_cw 0.99 0.92 1.06 1.001 8265.48 2975.70
is_A:pressure_partner_cw 0.98 0.90 1.06 1.001 7960.95 3211.18
is_B:pressure_partner_cw 0.96 0.90 1.03 1.001 7788.09 2798.46
is_A:day 0.89* 0.79 0.99 0.999 6857.79 3004.78
is_B:day 1.01 0.91 1.13 1.002 7204.53 3132.25
is_A:weartime_self_cw 1.00 1.00 1.00 1.001 3671.36 2458.77
is_B:weartime_self_cw 1.00 1.00 1.00 1.001 3917.54 2748.69
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.96 0.74 1.24 1.001 2660.19 2858.16
is_B:pressure_self_cb 0.57* 0.43 0.75 1.003 2662.74 3120.95
is_A:pressure_partner_cb 1.21 0.87 1.68 1.001 2634.87 3109.90
is_B:pressure_partner_cb 1.55* 1.25 1.92 1.002 2684.09 3001.05
is_A:weartime_self_cb 1.00 1.00 1.00 1.000 4690.54 3404.46
is_B:weartime_self_cb 1.00* 1.00 1.00 1.000 4172.07 3650.51
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.30 0.22 0.42 1.01 667.67 1367.32
sd(is_B) 0.39 0.29 0.52 1.01 519.76 1115.01
sd(is_A:pressure_self_cw) 0.05 0.00 0.15 1.00 2864.03 2463.77
sd(is_B:pressure_self_cw) 0.08 0.00 0.24 1.00 1698.22 1980.12
sd(is_A:pressure_partner_cw) 0.13 0.00 0.37 1.00 1112.61 1897.56
sd(is_B:pressure_partner_cw) 0.04 0.00 0.13 1.00 3045.19 1813.05
Additional Parameters
ar[1] 0.21 0.18 0.25 1.00 5907.45 2933.90
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.55 0.53 0.56 1.00 6026.96 2831.33

Compare Models

loo_compare(loo_ind, loo_apim)
##                          elpd_diff se_diff
## pressure_pa_obj_log_apim   0.0       0.0  
## pressure_pa_obj_log_ind  -72.7      14.3
bayes_R2(pressure_pa_obj_log_apim)
##    Estimate   Est.Error     Q2.5     Q97.5
## R2 0.183356 0.006415822 0.170636 0.1958512
bayes_R2(pressure_pa_obj_log_ind)
##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.1622604 0.006731342 0.1487327 0.175401

Affect ~ pressure

Indistinguishable

formula <- bf(
  aff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_aff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_aff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4821.6  63.5
## p_loo        58.1   3.3
## looic      9643.3 126.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_aff_ind, ask = FALSE)

summarize_brms(
  pressure_aff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.89* 4.67 5.10 1.001 745.62 1272.54
Within-Person Effects
pressure_self_cw -0.02 -0.10 0.05 1.000 5382.92 2754.52
pressure_partner_cw -0.01 -0.08 0.06 1.000 6190.02 2830.31
day -0.12 -0.27 0.03 1.002 5330.47 2728.44
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.12 -0.16 0.41 1.000 3197.40 2640.58
pressure_partner_cb 0.03 -0.24 0.32 1.001 3407.75 2855.72
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.47 0.78 1.01 788.16 1635.25
sd(pressure_self_cw) 0.10 0.00 0.26 1.00 832.42 1247.63
sd(pressure_partner_cw) 0.15 0.01 0.32 1.00 1122.34 1808.10
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 5392.27 2816.35
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 5374.82 2743.43

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_aff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_aff_apim)
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pressure_aff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4694.2  65.3
## p_loo       101.0   4.4
## looic      9388.3 130.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  238     
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_aff_apim, ask = FALSE)

summarize_brms(
  pressure_aff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 4.87* 4.58 5.13 1.002 618.64 1002.47
is_B 4.88* 4.60 5.17 1.003 605.34 891.66
Within-Person Effects
is_A:pressure_self_cw -0.03 -0.12 0.08 1.000 7705.26 2999.35
is_B:pressure_self_cw -0.05 -0.15 0.06 1.001 6987.03 3275.07
is_A:pressure_partner_cw -0.03 -0.13 0.06 1.000 6252.92 2711.34
is_B:pressure_partner_cw 0.01 -0.09 0.11 1.002 6068.39 2747.39
is_A:day -0.08 -0.25 0.10 1.004 6759.69 2583.51
is_B:day -0.13 -0.31 0.05 1.002 6070.96 2812.16
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.07 -0.34 0.47 1.000 2497.18 2873.84
is_B:pressure_self_cb 0.12 -0.33 0.56 1.001 2743.62 2687.95
is_A:pressure_partner_cb -0.02 -0.51 0.48 1.000 2586.37 2966.02
is_B:pressure_partner_cb -0.12 -0.48 0.24 1.000 2636.66 3049.75
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.79 0.62 1.01 1.00 903.34 2049.14
sd(is_B) 0.70 0.55 0.90 1.00 1001.21 1833.24
sd(is_A:pressure_self_cw) 0.10 0.00 0.31 1.00 1343.70 1872.04
sd(is_B:pressure_self_cw) 0.13 0.01 0.39 1.00 1767.18 2080.14
sd(is_A:pressure_partner_cw) 0.25 0.02 0.60 1.00 1117.76 1683.55
sd(is_B:pressure_partner_cw) 0.09 0.00 0.25 1.00 2168.54 1873.17
Additional Parameters
ar[1] 0.34 0.30 0.37 1.00 6535.64 2752.72
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.84 0.82 0.86 1.00 6234.22 2867.03

Compare Models

loo_compare(loo_ind, loo_apim)
##                   elpd_diff se_diff
## pressure_aff_apim    0.0       0.0 
## pressure_aff_ind  -127.5      18.6
bayes_R2(pressure_aff_apim)
##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.2349984 0.004573314 0.2260495 0.243759
bayes_R2(pressure_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2257293 0.004871448 0.2162224 0.2350853

Anticipatory Affect ~ pressure

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_anticipAff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7525.4  56.6
## p_loo        60.7   2.8
## looic     15050.8 113.2
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pressure_anticipAff_ind, ask = FALSE)

summarize_brms(
  pressure_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.46* 1.02 1.91 1.005 603.42 1144.57
Within-Person Effects
pressure_self_cw 0.04 -0.11 0.20 1.003 6227.17 2541.02
pressure_partner_cw 0.01 -0.14 0.16 1.000 6249.71 2958.94
day 0.01 -0.29 0.30 1.001 6655.31 2963.37
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.44 -0.10 0.96 1.001 2885.48 2768.50
pressure_partner_cb -0.32 -0.85 0.23 1.000 3041.60 2824.84
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.27 0.99 1.61 1.01 836.58 1818.51
sd(pressure_self_cw) 0.31 0.12 0.56 1.00 2121.08 2136.77
sd(pressure_partner_cw) 0.30 0.10 0.53 1.00 2000.95 1516.93
Additional Parameters
ar[1] 0.42 0.39 0.45 1.00 5240.60 2632.30
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.80 1.76 1.84 1.00 6329.26 3065.50

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_anticipAff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_anticipAff_apim)
## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pressure_anticipAff_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7424.0  56.4
## p_loo       106.8   4.5
## looic     14847.9 112.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3734  99.9%   390     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_anticipAff_apim, ask = FALSE)

summarize_brms(
  pressure_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.04* 0.43 1.62 1.005 821.47 1628.44
is_B 1.84* 1.26 2.44 1.006 689.95 1212.41
Within-Person Effects
is_A:pressure_self_cw -0.04 -0.25 0.17 1.001 8818.85 2730.05
is_B:pressure_self_cw 0.16 -0.06 0.38 1.000 8690.77 2508.23
is_A:pressure_partner_cw -0.10 -0.32 0.12 1.000 7399.40 3049.02
is_B:pressure_partner_cw 0.10 -0.10 0.31 1.003 8972.74 2544.64
is_A:day 0.22 -0.15 0.59 1.000 7710.82 3259.29
is_B:day -0.17 -0.56 0.21 1.002 6770.78 2914.78
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.11 -0.74 0.99 1.002 3186.28 2386.68
is_B:pressure_self_cb -0.56 -1.49 0.40 1.001 3634.63 2900.67
is_A:pressure_partner_cb 0.03 -1.08 1.11 1.003 3222.52 2956.47
is_B:pressure_partner_cb 0.44 -0.32 1.20 1.000 3554.36 3010.47
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 1.47 1.15 1.90 1.00 946.11 1892.12
sd(is_B) 1.49 1.16 1.92 1.01 1025.47 1516.31
sd(is_A:pressure_self_cw) 0.41 0.08 0.78 1.00 1295.17 1095.07
sd(is_B:pressure_self_cw) 0.33 0.05 0.71 1.00 1594.05 1563.60
sd(is_A:pressure_partner_cw) 0.41 0.13 0.81 1.00 3051.13 2709.20
sd(is_B:pressure_partner_cw) 0.24 0.01 0.63 1.00 1319.27 1629.73
Additional Parameters
ar[1] 0.34 0.30 0.37 1.00 6037.29 2884.95
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.74 1.70 1.78 1.01 8425.43 2231.59

Compare Models

loo_compare(loo_ind, loo_apim)
##                          elpd_diff se_diff
## pressure_anticipAff_apim    0.0       0.0 
## pressure_anticipAff_ind  -101.4      16.0
bayes_R2(pressure_anticipAff_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2202662 0.00455434 0.2110166 0.2288988
bayes_R2(pressure_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2112457 0.004909094 0.2016515 0.2206921

reactance ~ pressure

Indistinguishable

formula <- bf(
  reactance ~ 
    pressure_self_cw + pressure_partner_cw +
    pressure_self_cb + pressure_partner_cb +
    day +
    
    # Random effects
    (pressure_self_cw + pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_reactance_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pressure_reactance_ind)
## Warning: Found 9 observations with a pareto_k > 0.7 in model 'pressure_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1067.3 33.5
## p_loo        60.8  6.2
## looic      2134.7 67.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     747   98.8%   254     
##    (0.7, 1]   (bad)        8    1.1%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_reactance_ind, ask = FALSE)

summarize_brms(
  pressure_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pressure_ind,
  model_rows_random = model_rows_random_pressure_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.56* 0.33 0.79 1.001 1102.00 2248.74
Within-Person Effects
pressure_self_cw -0.02 -0.11 0.07 1.002 5289.22 3244.13
pressure_partner_cw -0.06 -0.18 0.06 1.002 6071.65 3117.34
day -0.13 -0.38 0.12 1.000 5474.68 3199.35
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 0.10 -0.20 0.39 1.001 2425.72 2868.45
pressure_partner_cb -0.05 -0.36 0.26 1.001 2329.16 2722.78
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.55 0.42 0.71 1.00 1575.28 2521.12
sd(pressure_self_cw) 0.44 0.27 0.66 1.00 1718.00 2543.78
sd(pressure_partner_cw) 0.17 0.01 0.51 1.00 964.20 1628.90
Additional Parameters
ar[1] 0.02 -0.06 0.10 1.00 4689.14 3011.29
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.95 0.90 1.00 1.00 4366.53 3027.93

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw +
    
    is_A:pressure_self_cb + is_B:pressure_self_cb + 
    is_A:pressure_partner_cb + is_B:pressure_partner_cb +
    
    is_A:day + is_B:day + 
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pressure_self_cw + is_B:pressure_self_cw + 
    is_A:pressure_partner_cw + is_B:pressure_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pressure_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  control = list(adapt_delta = 0.90, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pressure_reactance_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pressure_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pressure_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pressure_reactance_apim)
## Warning: Found 13 observations with a pareto_k > 0.7 in model 'pressure_reactance_apim'. We recommend to
## set 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1055.8 33.8
## p_loo        79.3  7.5
## looic      2111.5 67.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     743   98.3%   155     
##    (0.7, 1]   (bad)       13    1.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pressure_reactance_apim, ask = FALSE)

summarize_brms(
  pressure_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pressure_apim,
  model_rows_random = model_rows_random_pressure_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.29* 0.00 0.57 1.002 2749.60 3190.70
is_B 0.77* 0.46 1.07 1.001 2297.08 2841.29
Within-Person Effects
is_A:pressure_self_cw 0.09 -0.04 0.22 1.001 8356.25 2795.65
is_B:pressure_self_cw -0.12 -0.24 0.01 1.004 9228.00 3480.37
is_A:pressure_partner_cw 0.03 -0.16 0.23 1.003 7949.62 3269.05
is_B:pressure_partner_cw -0.10 -0.25 0.04 1.003 10265.45 3046.31
is_A:day 0.12 -0.25 0.49 1.004 6959.77 3167.17
is_B:day -0.25 -0.63 0.15 1.000 6353.28 3401.10
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 0.51 -0.03 1.05 1.001 3906.50 3115.36
is_B:pressure_self_cb -0.33 -0.96 0.28 1.001 3785.08 3021.86
is_A:pressure_partner_cb -0.69 -1.40 0.01 1.000 3981.60 3159.78
is_B:pressure_partner_cb 0.33 -0.14 0.83 1.000 4266.32 3111.87
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA
Random Effects
sd(is_A) 0.55 0.40 0.74 1.00 2046.58 3141.28
sd(is_B) 0.50 0.35 0.69 1.00 2747.64 3247.90
sd(is_A:pressure_self_cw) 0.42 0.23 0.70 1.00 2428.64 3141.31
sd(is_B:pressure_self_cw) 0.48 0.25 0.84 1.00 2588.51 3414.47
sd(is_A:pressure_partner_cw) 0.36 0.03 0.85 1.00 1720.32 1800.65
sd(is_B:pressure_partner_cw) 0.10 0.00 0.35 1.00 2722.46 2590.97
Additional Parameters
ar[1] -0.02 -0.10 0.07 1.00 5172.64 3003.46
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.92 0.87 0.97 1.00 4939.79 3108.00

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pressure_reactance_apim   0.0       0.0  
## pressure_reactance_ind  -11.6       6.1
bayes_R2(pressure_reactance_apim)
##    Estimate   Est.Error     Q2.5     Q97.5
## R2 0.243469 0.009813226 0.223816 0.2621822
bayes_R2(pressure_reactance_ind)
##     Estimate   Est.Error     Q2.5     Q97.5
## R2 0.2322911 0.009985695 0.212067 0.2511857

Report models for pressure side by side

Indistinguishable Models

pressure_all_ind_models <- report_side_by_side(
  pressure_pa_sub_ind,
  pressure_pa_obj_log_ind,
  pressure_aff_ind,
  pressure_anticipAff_ind,
  pressure_reactance_ind,
  
  model_rows_random = model_rows_random_pressure_ind,
  model_rows_fixed = model_rows_fixed_pressure_ind
) 
## [1] "pressure_pa_sub_ind"
## [1] "pressure_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pressure_aff_ind"
## [1] "pressure_anticipAff_ind"
## [1] "pressure_reactance_ind"
# pretty printing
pressure_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pressure_pa_sub_ind 95% CI pressure_pa_sub_ind exp(Est.) pressure_pa_obj_log_ind 95% CI pressure_pa_obj_log_ind b pressure_aff_ind 95% CI pressure_aff_ind b pressure_anticipAff_ind 95% CI pressure_anticipAff_ind b pressure_reactance_ind 95% CI pressure_reactance_ind
Intercept 24.44* [18.79, 32.14] 118.46* [106.30, 131.66] 4.89* [ 4.67, 5.10] 1.46* [ 1.02, 1.91] 0.56* [ 0.33, 0.79]
Within-Person Effects
pressure_self_cw 0.98 [ 0.79, 1.23] 1.00 [ 0.95, 1.05] -0.02 [-0.10, 0.05] 0.04 [-0.11, 0.20] -0.02 [-0.11, 0.07]
pressure_partner_cw 0.86 [ 0.69, 1.10] 0.96 [ 0.91, 1.02] -0.01 [-0.08, 0.06] 0.01 [-0.14, 0.16] -0.06 [-0.18, 0.06]
day 1.18 [ 0.86, 1.62] 0.97 [ 0.89, 1.06] -0.12 [-0.27, 0.03] 0.01 [-0.29, 0.30] -0.13 [-0.38, 0.12]
weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
pressure_self_cb 1.09 [ 0.66, 1.81] 0.78* [ 0.67, 0.93] 0.12 [-0.16, 0.41] 0.44 [-0.10, 0.96] 0.10 [-0.20, 0.39]
pressure_partner_cb 0.89 [ 0.55, 1.47] 1.18* [ 1.02, 1.37] 0.03 [-0.24, 0.32] -0.32 [-0.85, 0.23] -0.05 [-0.36, 0.26]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.63 [ 0.47, 0.84] 0.29 [0.22, 0.37] 0.60 [0.47, 0.78] 1.27 [0.99, 1.61] 0.55 [ 0.42, 0.71]
sd(pressure_self_cw) 0.13 [ 0.01, 0.39] 0.05 [0.00, 0.14] 0.10 [0.00, 0.26] 0.31 [0.12, 0.56] 0.44 [ 0.27, 0.66]
sd(pressure_partner_cw) 0.16 [ 0.00, 0.47] 0.03 [0.00, 0.10] 0.15 [0.01, 0.32] 0.30 [0.10, 0.53] 0.17 [ 0.01, 0.51]
Additional Parameters
ar[1] 0.01 [-0.93, 0.94] 0.28 [0.25, 0.32] 0.45 [0.42, 0.48] 0.42 [0.39, 0.45] 0.02 [-0.06, 0.10]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.13 [ 0.12, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.04 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.56 [0.55, 0.58] 0.87 [0.85, 0.89] 1.80 [1.76, 1.84] 0.95 [ 0.90, 1.00]
# prepare export to excel
all_models$Pressure_ind <- prepare_for_export(pressure_all_ind_models)

APIMs

pressure_all_apim_models <- report_side_by_side(
  pressure_pa_sub_apim,
  pressure_pa_obj_log_apim,
  pressure_aff_apim,
  pressure_anticipAff_apim,
  pressure_reactance_apim,
  
  model_rows_random = model_rows_random_pressure_apim,
  model_rows_fixed = model_rows_fixed_pressure_apim
) 
## [1] "pressure_pa_sub_apim"
## [1] "pressure_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pressure_aff_apim"
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## [1] "pressure_anticipAff_apim"
## [1] "pressure_reactance_apim"
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
# pretty printing
pressure_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pressure_pa_sub_apim 95% CI pressure_pa_sub_apim exp(Est.) pressure_pa_obj_log_apim 95% CI pressure_pa_obj_log_apim b pressure_aff_apim 95% CI pressure_aff_apim b pressure_anticipAff_apim 95% CI pressure_anticipAff_apim b pressure_reactance_apim 95% CI pressure_reactance_apim
is_A 29.61* [20.61, 41.74] 127.96* [106.77, 150.24] 4.87* [ 4.58, 5.13] 1.04* [ 0.43, 1.62] 0.29* [ 0.00, 0.57]
is_B 19.98* [14.03, 28.42] 123.39* [104.27, 145.22] 4.88* [ 4.60, 5.17] 1.84* [ 1.26, 2.44] 0.77* [ 0.46, 1.07]
Within-Person Effects
is_A:pressure_self_cw 0.78 [ 0.59, 1.07] 1.01 [ 0.94, 1.09] -0.03 [-0.12, 0.08] -0.04 [-0.25, 0.17] 0.09 [-0.04, 0.22]
is_B:pressure_self_cw 1.26 [ 0.91, 1.88] 0.99 [ 0.92, 1.06] -0.05 [-0.15, 0.06] 0.16 [-0.06, 0.38] -0.12 [-0.24, 0.01]
is_A:pressure_partner_cw 1.02 [ 0.76, 1.45] 0.98 [ 0.90, 1.06] -0.03 [-0.13, 0.06] -0.10 [-0.32, 0.12] 0.03 [-0.16, 0.23]
is_B:pressure_partner_cw 0.63* [ 0.44, 0.95] 0.96 [ 0.90, 1.03] 0.01 [-0.09, 0.11] 0.10 [-0.10, 0.31] -0.10 [-0.25, 0.04]
is_A:day 1.07 [ 0.69, 1.68] 0.89* [ 0.79, 0.99] -0.08 [-0.25, 0.10] 0.22 [-0.15, 0.59] 0.12 [-0.25, 0.49]
is_B:day 1.20 [ 0.78, 1.81] 1.01 [ 0.91, 1.13] -0.13 [-0.31, 0.05] -0.17 [-0.56, 0.21] -0.25 [-0.63, 0.15]
is_A:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cw NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
is_A:pressure_self_cb 1.08 [ 0.49, 2.36] 0.96 [ 0.74, 1.24] 0.07 [-0.34, 0.47] 0.11 [-0.74, 0.99] 0.51 [-0.03, 1.05]
is_B:pressure_self_cb 1.16 [ 0.45, 3.08] 0.57* [ 0.43, 0.75] 0.12 [-0.33, 0.56] -0.56 [-1.49, 0.40] -0.33 [-0.96, 0.28]
is_A:pressure_partner_cb 0.79 [ 0.29, 2.13] 1.21 [ 0.87, 1.68] -0.02 [-0.51, 0.48] 0.03 [-1.08, 1.11] -0.69 [-1.40, 0.01]
is_B:pressure_partner_cb 0.84 [ 0.41, 1.76] 1.55* [ 1.25, 1.92] -0.12 [-0.48, 0.24] 0.44 [-0.32, 1.20] 0.33 [-0.14, 0.83]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
is_B:barriers_self_cb NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(is_A) 0.66 [ 0.47, 0.90] 0.30 [0.22, 0.42] 0.79 [0.62, 1.01] 1.47 [1.15, 1.90] 0.55 [ 0.40, 0.74]
sd(is_B) 0.57 [ 0.37, 0.80] 0.39 [0.29, 0.52] 0.70 [0.55, 0.90] 1.49 [1.16, 1.92] 0.50 [ 0.35, 0.69]
sd(is_A:pressure_self_cw) 0.21 [ 0.01, 0.70] 0.05 [0.00, 0.15] 0.10 [0.00, 0.31] 0.41 [0.08, 0.78] 0.42 [ 0.23, 0.70]
sd(is_B:pressure_self_cw) 0.21 [ 0.01, 0.72] 0.08 [0.00, 0.24] 0.13 [0.01, 0.39] 0.33 [0.05, 0.71] 0.48 [ 0.25, 0.84]
sd(is_A:pressure_partner_cw) 0.25 [ 0.01, 0.85] 0.13 [0.00, 0.37] 0.25 [0.02, 0.60] 0.41 [0.13, 0.81] 0.36 [ 0.03, 0.85]
sd(is_B:pressure_partner_cw) 0.24 [ 0.01, 0.82] 0.04 [0.00, 0.13] 0.09 [0.00, 0.25] 0.24 [0.01, 0.63] 0.10 [ 0.00, 0.35]
Additional Parameters
ar[1] 0.02 [-0.94, 0.94] 0.21 [0.18, 0.25] 0.34 [0.30, 0.37] 0.34 [0.30, 0.37] -0.02 [-0.10, 0.07]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.13 [ 0.13, 0.14] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.04 [ 0.00, 0.12] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.55 [0.53, 0.56] 0.84 [0.82, 0.86] 1.74 [1.70, 1.78] 0.92 [ 0.87, 0.97]
# prepare export to excel
all_models$Pressure_apim <- prepare_for_export(pressure_all_apim_models)

Modelling Pushing

# For indistinguishable Dyads
model_rows_fixed_pushing_ind <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw', 
    'barriers_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    'barriers_self_cb'
  )
  
model_rows_random_pushing_ind <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)


# For APIMs

model_rows_fixed_pushing_apim <- c(
    'is_A',
    'is_B',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'is_A:pushing_self_cw',
    'is_B:pushing_self_cw', 
    'is_A:pushing_partner_cw', 
    'is_B:pushing_partner_cw',
    'is_A:day', 
    'is_B:day', 
    'is_A:weartime_self_cw', 
    'is_B:weartime_self_cw', 
    'is_A:barriers_self_cw',
    'is_B:barriers_self_cw',
    
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'is_A:pushing_self_cb',
    'is_B:pushing_self_cb', 
    'is_A:pushing_partner_cb', 
    'is_B:pushing_partner_cb',
    'is_A:weartime_self_cb', 
    'is_B:weartime_self_cb', 
    'is_A:barriers_self_cb',
    'is_B:barriers_self_cb'
  )
  

model_rows_random_pushing_apim <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(is_A)',
  'sd(is_B)',
  'sd(is_A:pushing_self_cw)',
  'sd(is_B:pushing_self_cw)', 
  'sd(is_A:pushing_partner_cw)', 
  'sd(is_B:pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'nu',
  'shape',
  'sderr',
  'sigma'
)

Subjective MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_sub_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_sub_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_sub_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_sub_ind)
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pushing_pa_sub_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6316.5  67.5
## p_loo        33.5   2.8
## looic     12633.0 135.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1339  99.9%   242     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_sub_ind, ask = FALSE)

summarize_brms(
  pushing_pa_sub_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 49.28* 37.64 64.23 1.001 1655.14 2441.46
Within-Person Effects
pushing_self_cw 1.00 0.90 1.11 1.000 5361.40 3118.25
pushing_partner_cw 1.02 0.92 1.13 1.000 5970.95 2971.70
day 0.86 0.64 1.15 1.001 6379.73 2842.66
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 1.02 0.95 1.08 1.000 7012.08 3036.44
Between-Person Effects
pushing_self_cb 0.59 0.34 1.03 1.000 3443.47 3192.15
pushing_partner_cb 1.50 0.91 2.49 1.000 5389.19 3197.07
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 1.03 0.90 1.17 1.001 3903.33 2745.51
Random Effects
sd(Intercept) 0.63 0.44 0.87 1.00 1301.04 2033.20
sd(pushing_self_cw) 0.31 0.13 0.51 1.00 1453.47 1985.91
sd(pushing_partner_cw) 0.28 0.11 0.48 1.00 1436.65 1793.03
Additional Parameters
ar[1] 0.04 -0.93 0.95 1.00 5529.16 2688.10
nu NA NA NA NA NA NA
shape 0.43 0.40 0.47 1.00 5112.46 3066.21
sderr 0.04 0.00 0.12 1.00 2385.84 2277.18
sigma NA NA NA NA NA NA

APIM

formula <- bf(
  pa_sub ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +
    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 20)", class = "shape"), 
  brms::set_prior("cauchy(0, 10)", class='sderr')
)

df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_sub_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = brms::negbinomial(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_sub_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_sub_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_sub_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_sub_apim)
## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_pa_sub_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6323.5  67.5
## p_loo        43.1   2.8
## looic     12647.1 135.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1338  99.9%   292     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_sub_apim, ask = FALSE)

summarize_brms(
  pushing_pa_sub_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
IRR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 64.43* 46.15 89.96 1.000 2298.52 2733.60
is_B 35.31* 25.01 49.86 1.000 2570.13 3174.92
Within-Person Effects
is_A:pushing_self_cw 1.05 0.92 1.21 1.001 7211.92 2859.27
is_B:pushing_self_cw 0.92 0.79 1.09 1.002 6161.35 3421.44
is_A:pushing_partner_cw 1.03 0.89 1.20 1.000 7086.77 3063.78
is_B:pushing_partner_cw 1.00 0.87 1.15 1.001 6245.30 3188.98
is_A:day 0.73 0.48 1.12 1.000 5634.86 3152.56
is_B:day 1.07 0.69 1.65 1.000 5717.33 3215.41
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 1.03 0.93 1.14 1.002 6587.34 2507.94
is_B:barriers_self_cw 0.99 0.91 1.08 1.000 7715.28 3127.70
Between-Person Effects
is_A:pushing_self_cb 0.47 0.19 1.15 1.000 3547.21 3588.53
is_B:pushing_self_cb 0.73 0.35 1.52 1.000 5000.77 3408.97
is_A:pushing_partner_cb 1.31 0.68 2.60 1.001 5162.93 3354.12
is_B:pushing_partner_cb 2.03 0.84 4.97 1.001 4122.35 3312.50
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 1.08 0.89 1.31 1.000 4317.48 3377.06
is_B:barriers_self_cb 1.04 0.82 1.30 1.001 3580.54 2922.44
Random Effects
sd(is_A) 0.63 0.39 0.97 1.00 1165.20 1894.07
sd(is_B) 0.47 0.29 0.68 1.00 1978.49 2586.86
sd(is_A:pushing_self_cw) 0.12 0.01 0.33 1.00 1640.62 2109.66
sd(is_B:pushing_self_cw) 0.46 0.14 0.85 1.00 1295.91 1465.64
sd(is_A:pushing_partner_cw) 0.43 0.11 0.84 1.00 1172.68 1373.02
sd(is_B:pushing_partner_cw) 0.11 0.01 0.30 1.00 2517.62 2903.67
Additional Parameters
ar[1] 0.04 -0.92 0.94 1.01 5063.36 2243.66
nu NA NA NA NA NA NA
shape 0.43 0.40 0.47 1.00 4753.53 2795.10
sderr 0.05 0.00 0.13 1.00 2691.83 2478.37
sigma NA NA NA NA NA NA

Compare models

loo_compare(loo_ind, loo_apim)
##                     elpd_diff se_diff
## pushing_pa_sub_ind   0.0       0.0   
## pushing_pa_sub_apim -7.0       3.6
bayes_R2(pushing_pa_sub_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2654618 0.07201488 0.1666211 0.4672379
bayes_R2(pushing_pa_sub_ind)
##    Estimate  Est.Error      Q2.5     Q97.5
## R2 0.213283 0.04997866 0.1402257 0.3408443

Device Based MVPA ~ pushing

Indistinguishable

formula <- bf(
  pa_obj_log ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + weartime_self_cw + weartime_self_cb +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_obj_log_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_obj_log_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_obj_log_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_pa_obj_log_ind)
print(loo_ind)
## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -981.7 34.0
## p_loo        64.5  4.7
## looic      1963.4 68.0
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_obj_log_ind, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = T) %>%
  print_df() %>%
  packing_ind()
## Warning in summarize_brms(pushing_pa_obj_log_ind, model_rows_fixed = model_rows_fixed_pushing_ind, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 127.21* 111.99 144.82 1.000 1752.75 2358.30
Within-Person Effects
pushing_self_cw 1.00 0.97 1.04 1.001 7809.82 2778.98
pushing_partner_cw 0.99 0.96 1.03 1.000 7975.04 2950.37
day 1.00 0.88 1.14 1.000 7752.54 3200.40
weartime_self_cw 1.00 1.00 1.00 1.001 3820.73 2574.65
barriers_self_cw 0.98 0.96 1.01 1.001 7182.40 2846.49
Between-Person Effects
pushing_self_cb 0.82 0.64 1.07 1.002 4513.32 3468.94
pushing_partner_cb 0.86 0.66 1.12 1.001 4752.33 3505.51
weartime_self_cb 1.00 1.00 1.00 1.002 5376.81 3249.84
barriers_self_cb 0.99 0.93 1.05 1.002 4331.99 3156.88
Random Effects
sd(Intercept) 0.30 0.22 0.41 1.00 1217.79 2263.40
sd(pushing_self_cw) 0.10 0.04 0.17 1.00 1221.10 1899.51
sd(pushing_partner_cw) 0.04 0.00 0.11 1.00 1322.25 1982.92
Additional Parameters
ar[1] 0.27 0.21 0.33 1.00 5406.99 2539.64
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.53 0.51 0.55 1.00 6367.72 2955.61

apim

formula <- bf(
  pa_obj_log ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:weartime_self_cw + is_B:weartime_self_cw + 
    is_A:weartime_self_cb + is_B:weartime_self_cb +
    
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b"),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_A'),
  brms::set_prior("normal(0, 50)", class = "b", coef = 'is_B'),
  
  brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
  
  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_pa_obj_log_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.99),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_pa_obj_log_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_pa_obj_log_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_pa_obj_log_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_pa_obj_log_apim)
## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_pa_obj_log_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 1203 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -972.6 35.2
## p_loo       106.4  7.1
## looic      1945.3 70.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1201  99.8%   86      
##    (0.7, 1]   (bad)         2   0.2%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_pa_obj_log_apim, ask = FALSE)

summarize_brms(
  pushing_pa_obj_log_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = T) %>%
  print_df() %>%
  packing_apim()
## Warning in summarize_brms(pushing_pa_obj_log_apim, model_rows_fixed = model_rows_fixed_pushing_apim, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 129.10* 106.65 155.78 1.003 831.15 1018.63
is_B 126.93* 104.96 153.25 1.004 910.22 1647.17
Within-Person Effects
is_A:pushing_self_cw 0.99 0.95 1.04 1.001 5241.24 2970.02
is_B:pushing_self_cw 1.02 0.96 1.08 1.000 5683.20 2956.87
is_A:pushing_partner_cw 1.03 0.97 1.09 1.000 4955.72 2932.87
is_B:pushing_partner_cw 0.98 0.93 1.02 1.000 6065.71 3072.27
is_A:day 1.08 0.91 1.28 1.000 4253.08 3276.88
is_B:day 0.94 0.78 1.12 1.000 3761.56 2977.15
is_A:weartime_self_cw 1.00 1.00 1.00 1.001 4018.46 2772.97
is_B:weartime_self_cw 1.00 1.00 1.00 1.002 4229.94 2751.28
is_A:barriers_self_cw 1.01 0.98 1.04 1.000 5530.29 3127.49
is_B:barriers_self_cw 0.98 0.95 1.00 1.000 5543.34 2635.67
Between-Person Effects
is_A:pushing_self_cb 0.81 0.56 1.19 1.000 3048.72 2968.96
is_B:pushing_self_cb 0.77 0.52 1.13 1.001 2956.11 2972.07
is_A:pushing_partner_cb 1.12 0.77 1.64 1.001 3227.02 3033.76
is_B:pushing_partner_cb 0.74 0.48 1.16 1.002 2533.24 2648.69
is_A:weartime_self_cb 1.00 1.00 1.00 1.000 3410.36 3401.54
is_B:weartime_self_cb 1.00 1.00 1.00 1.000 4222.92 3290.03
is_A:barriers_self_cb 1.07 0.99 1.17 1.001 2854.71 2630.65
is_B:barriers_self_cb 0.98 0.89 1.09 1.001 2917.54 3217.70
Random Effects
sd(is_A) 0.34 0.24 0.47 1.00 1058.13 2038.29
sd(is_B) 0.34 0.23 0.48 1.00 857.30 1343.85
sd(is_A:pushing_self_cw) 0.09 0.01 0.18 1.01 911.24 911.67
sd(is_B:pushing_self_cw) 0.11 0.01 0.20 1.00 922.39 638.42
sd(is_A:pushing_partner_cw) 0.11 0.02 0.20 1.00 1053.84 885.89
sd(is_B:pushing_partner_cw) 0.03 0.00 0.09 1.00 2359.01 2068.39
Additional Parameters
ar[1] 0.21 0.14 0.28 1.00 3704.50 3396.22
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.52 0.50 0.54 1.00 3544.78 2792.99

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pushing_pa_obj_log_apim  0.0       0.0   
## pushing_pa_obj_log_ind  -9.1       7.6
bayes_R2(pushing_pa_obj_log_apim)
##    Estimate  Est.Error     Q2.5     Q97.5
## R2 0.210763 0.01165598 0.187023 0.2326488
bayes_R2(pushing_pa_obj_log_ind)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.1859197  0.012287 0.1615149 0.2094954

Affect ~ pushing

Indistinguishable

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_aff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_aff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_aff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_aff_ind)
print(loo_ind)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1769.7 41.4
## p_loo        55.0  3.4
## looic      3539.5 82.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_aff_ind, ask = FALSE)

summarize_brms(
  pushing_aff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 5.00* 4.77 5.24 1.008 705.32 1279.78
Within-Person Effects
pushing_self_cw -0.03 -0.08 0.03 1.004 5478.74 3016.82
pushing_partner_cw 0.03 -0.02 0.09 1.001 5545.96 2960.22
day -0.04 -0.24 0.15 1.001 4709.38 2754.74
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 0.00 -0.04 0.03 1.001 5371.78 3127.61
Between-Person Effects
pushing_self_cb -0.07 -0.49 0.33 1.000 2734.11 3086.61
pushing_partner_cb 0.22 -0.16 0.60 1.000 3505.84 2982.04
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb -0.08 -0.18 0.02 1.001 3245.62 2991.88
Random Effects
sd(Intercept) 0.62 0.47 0.81 1.00 971.33 1624.25
sd(pushing_self_cw) 0.07 0.01 0.15 1.00 1114.99 1379.33
sd(pushing_partner_cw) 0.12 0.03 0.22 1.00 1204.91 725.56
Additional Parameters
ar[1] 0.27 0.22 0.33 1.00 4812.62 3064.25
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.89 0.86 0.92 1.00 5261.68 2709.15

APIM

formula <- bf(
  aff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
  
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_aff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_aff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_aff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_aff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_aff_apim)
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1740.0 41.3
## p_loo        98.2  5.7
## looic      3480.0 82.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_aff_apim, ask = FALSE)

summarize_brms(
  pushing_aff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 5.00* 4.71 5.29 1.001 1540.42 2422.92
is_B 4.90* 4.59 5.21 1.001 1294.81 2299.76
Within-Person Effects
is_A:pushing_self_cw -0.01 -0.08 0.07 1.001 7720.22 2873.36
is_B:pushing_self_cw -0.06 -0.14 0.03 1.000 8447.37 3343.11
is_A:pushing_partner_cw 0.00 -0.08 0.09 1.000 7345.86 2894.63
is_B:pushing_partner_cw 0.05 -0.02 0.13 1.002 8307.39 3042.29
is_A:day -0.06 -0.33 0.21 1.000 6880.69 3113.49
is_B:day -0.01 -0.26 0.25 1.001 7737.84 3345.38
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw -0.02 -0.07 0.04 1.003 8714.70 2540.72
is_B:barriers_self_cw 0.02 -0.03 0.06 1.001 7904.22 2670.31
Between-Person Effects
is_A:pushing_self_cb 0.42 -0.22 1.08 1.002 3505.56 3356.73
is_B:pushing_self_cb -0.28 -0.84 0.26 1.000 4412.31 3472.92
is_A:pushing_partner_cb -0.07 -0.61 0.49 1.000 4266.48 3331.04
is_B:pushing_partner_cb 0.41 -0.23 1.06 1.000 4079.82 3420.17
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 0.01 -0.14 0.16 1.000 3759.93 3437.66
is_B:barriers_self_cb -0.07 -0.25 0.11 1.000 2812.54 2230.19
Random Effects
sd(is_A) 0.79 0.61 1.02 1.00 1664.25 2441.20
sd(is_B) 0.68 0.51 0.89 1.00 1527.28 2451.03
sd(is_A:pushing_self_cw) 0.07 0.00 0.18 1.00 2328.67 2045.82
sd(is_B:pushing_self_cw) 0.12 0.01 0.27 1.00 1568.23 1969.97
sd(is_A:pushing_partner_cw) 0.17 0.04 0.32 1.00 1934.22 1642.01
sd(is_B:pushing_partner_cw) 0.10 0.01 0.24 1.00 1728.04 2285.82
Additional Parameters
ar[1] 0.18 0.12 0.24 1.00 4212.94 2703.09
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 0.85 0.82 0.89 1.00 5806.11 2871.44

Compare Models

loo_compare(loo_ind, loo_apim)
##                  elpd_diff se_diff
## pushing_aff_apim   0.0       0.0  
## pushing_aff_ind  -29.7      10.3
bayes_R2(pushing_aff_apim)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2077469 0.007994099 0.1918707 0.2232979
bayes_R2(pushing_aff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.1898545 0.008693353 0.1725112 0.2062479

Anticipatory Affect ~ pushing

Indistinguishable

formula <- bf(
  anticipAff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day + 
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=-5, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_anticipAff_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_anticipAff_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_anticipAff_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_anticipAff_ind)
print(loo_ind)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2770.8 30.9
## p_loo        51.8  2.6
## looic      5541.5 61.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.6]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
plot(pushing_anticipAff_ind, ask = FALSE)

summarize_brms(
  pushing_anticipAff_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.91* 1.44 2.37 1.007 804.61 1509.41
Within-Person Effects
pushing_self_cw 0.01 -0.11 0.12 1.001 5355.07 2643.69
pushing_partner_cw -0.07 -0.19 0.05 1.003 4170.72 3032.32
day -0.25 -0.66 0.19 1.000 4744.55 3205.48
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw -0.01 -0.09 0.06 1.000 5164.31 2685.59
Between-Person Effects
pushing_self_cb -0.34 -1.22 0.55 1.001 2616.54 2962.99
pushing_partner_cb -0.05 -0.91 0.82 1.000 3466.69 2656.87
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 0.17 -0.06 0.39 1.002 2864.71 2851.86
Random Effects
sd(Intercept) 1.22 0.92 1.60 1.00 1047.94 2128.53
sd(pushing_self_cw) 0.15 0.01 0.34 1.00 1188.27 1230.81
sd(pushing_partner_cw) 0.17 0.01 0.39 1.00 1192.06 1434.02
Additional Parameters
ar[1] 0.34 0.28 0.39 1.00 4137.46 3139.97
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.88 1.80 1.95 1.00 4819.75 2935.38

APIM

formula <- bf(
  anticipAff ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_anticipAff_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_anticipAff_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_anticipAff_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_anticipAff_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_anticipAff_apim)
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'pushing_anticipAff_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 1340 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2738.5 31.8
## p_loo        89.1  4.4
## looic      5477.0 63.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     1339  99.9%   737     
##    (0.7, 1]   (bad)         1   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_anticipAff_apim, ask = FALSE)

summarize_brms(
  pushing_anticipAff_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 1.60* 0.96 2.23 1.000 2336.15 2468.72
is_B 2.02* 1.41 2.65 1.002 2400.99 2975.07
Within-Person Effects
is_A:pushing_self_cw 0.04 -0.11 0.19 1.000 8718.69 3060.17
is_B:pushing_self_cw -0.07 -0.25 0.10 1.001 8885.25 3221.91
is_A:pushing_partner_cw 0.00 -0.18 0.17 1.001 9114.76 2944.38
is_B:pushing_partner_cw -0.09 -0.25 0.06 1.001 7658.91 3226.08
is_A:day 0.06 -0.51 0.64 1.001 7604.15 3302.30
is_B:day -0.55 -1.13 0.05 1.000 9013.42 2681.08
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 0.01 -0.10 0.12 1.000 10377.14 3052.80
is_B:barriers_self_cw -0.03 -0.13 0.06 1.001 8435.68 2823.84
Between-Person Effects
is_A:pushing_self_cb 1.01 -0.47 2.54 1.000 4815.92 3066.66
is_B:pushing_self_cb -0.80 -1.96 0.41 1.002 5833.57 3178.21
is_A:pushing_partner_cb -1.21 -2.51 0.05 1.000 5793.34 3572.47
is_B:pushing_partner_cb 1.16 -0.36 2.58 1.001 5095.68 3114.79
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 0.21 -0.11 0.54 1.000 5009.92 3312.75
is_B:barriers_self_cb 0.04 -0.35 0.42 1.001 4198.81 3167.66
Random Effects
sd(is_A) 1.43 1.08 1.88 1.00 1466.89 2692.09
sd(is_B) 1.52 1.15 1.98 1.00 2066.01 2792.50
sd(is_A:pushing_self_cw) 0.15 0.01 0.37 1.00 2017.59 2134.81
sd(is_B:pushing_self_cw) 0.19 0.01 0.49 1.00 1942.39 2438.13
sd(is_A:pushing_partner_cw) 0.28 0.02 0.60 1.00 1525.13 1808.11
sd(is_B:pushing_partner_cw) 0.11 0.00 0.31 1.00 2727.86 2989.60
Additional Parameters
ar[1] 0.24 0.18 0.30 1.00 5413.34 2912.22
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.80 1.73 1.88 1.00 7463.90 2647.65

Compare Models

loo_compare(loo_ind, loo_apim)
##                         elpd_diff se_diff
## pushing_anticipAff_apim   0.0       0.0  
## pushing_anticipAff_ind  -32.3       9.7
bayes_R2(pushing_anticipAff_apim)
##     Estimate   Est.Error      Q2.5    Q97.5
## R2 0.2212212 0.008630534 0.2035636 0.237974
bayes_R2(pushing_anticipAff_ind)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.2016732 0.009264627 0.1836161 0.2195859

reactance ~ pushing

Indistinguishable

formula <- bf(
  reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    day +
    barriers_self_cw + barriers_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
  
  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_ind <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_reactance_ind")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_reactance_ind, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_reactance_ind)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_ind <- loo(pushing_reactance_ind)
## Warning: Found 2 observations with a pareto_k > 0.7 in model 'pushing_reactance_ind'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_ind)
## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -597.3 24.5
## p_loo        38.4  4.8
## looic      1194.6 49.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.0]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     401   99.5%   301     
##    (0.7, 1]   (bad)        2    0.5%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_reactance_ind, ask = FALSE)

summarize_brms(
  pushing_reactance_ind, 
  model_rows_fixed = model_rows_fixed_pushing_ind,
  model_rows_random = model_rows_random_pushing_ind,
  exponentiate = F) %>%
  print_df() %>%
  packing_ind()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.48* 0.25 0.71 1.002 6101.16 3321.85
Within-Person Effects
pushing_self_cw -0.02 -0.10 0.07 1.002 7139.75 3154.64
pushing_partner_cw 0.02 -0.07 0.11 1.001 8174.34 3071.49
day 0.27 -0.12 0.66 1.003 8161.51 2643.71
weartime_self_cw NA NA NA NA NA NA
barriers_self_cw 0.05 -0.03 0.12 1.003 8452.29 3346.46
Between-Person Effects
pushing_self_cb -0.71* -1.22 -0.20 1.000 5670.61 3274.64
pushing_partner_cb 0.27 -0.31 0.87 1.000 5673.75 3392.01
weartime_self_cb NA NA NA NA NA NA
barriers_self_cb 0.01 -0.17 0.18 1.002 4772.32 3303.19
Random Effects
sd(Intercept) 0.17 0.01 0.38 1.00 1653.80 2222.19
sd(pushing_self_cw) 0.23 0.11 0.37 1.00 1449.65 1987.04
sd(pushing_partner_cw) 0.08 0.00 0.23 1.00 1770.45 1884.97
Additional Parameters
ar[1] 0.16 0.03 0.28 1.00 5132.70 3219.90
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.01 0.94 1.09 1.00 5895.25 3394.64

APIM

formula <- bf(
  reactance ~ 0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw +
    
    is_A:pushing_self_cb + is_B:pushing_self_cb + 
    is_A:pushing_partner_cb + is_B:pushing_partner_cb +
    
    is_A:day + is_B:day + 
    is_A:barriers_self_cw + is_B:barriers_self_cw + 
    is_A:barriers_self_cb + is_B:barriers_self_cb +

    
    # Random effects
    (0 + 
    is_A + is_B + # Intercepts
    is_A:pushing_self_cw + is_B:pushing_self_cw + 
    is_A:pushing_partner_cw + is_B:pushing_partner_cw | coupleID),

  autocor = ~ ar(time = day, gr = coupleID:AorB, p = 1)
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b"),
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_A'), 
  brms::set_prior("normal(0, 20)", class = "b", coef = 'is_B'), 

  brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),

  brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
  brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)


df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pushing_reactance_apim <- brm(
  formula, 
  prior = prior1,
  data = df_minimal, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  #iter = 8000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", "as_preregistered_pushing_reactance_apim")
)
## Warning: Rows containing NAs were excluded from the model.
pp_check(pushing_reactance_apim, type='hist')
## Using 10 posterior draws for ppc type 'hist' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(pushing_reactance_apim)
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_apim <- loo(pushing_reactance_apim)
## Warning: Found 4 observations with a pareto_k > 0.7 in model 'pushing_reactance_apim'. We recommend to set
## 'moment_match = TRUE' in order to perform moment matching for problematic observations.
print(loo_apim)
## 
## Computed from 4000 by 403 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -603.3 24.4
## p_loo        58.7  6.8
## looic      1206.6 48.7
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     399   99.0%   118     
##    (0.7, 1]   (bad)        4    1.0%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
plot(pushing_reactance_apim, ask = FALSE)

summarize_brms(
  pushing_reactance_apim, 
  model_rows_fixed = model_rows_fixed_pushing_apim,
  model_rows_random = model_rows_random_pushing_apim,
  exponentiate = F) %>%
  print_df() %>%
  packing_apim()
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
is_A 0.30 -0.01 0.60 1.001 3000.47 2975.13
is_B 0.74* 0.41 1.08 1.004 2659.02 2802.83
Within-Person Effects
is_A:pushing_self_cw 0.00 -0.12 0.11 1.002 4557.74 3278.77
is_B:pushing_self_cw -0.05 -0.18 0.08 1.000 4887.79 3100.72
is_A:pushing_partner_cw 0.00 -0.13 0.13 1.001 5093.91 3230.72
is_B:pushing_partner_cw 0.05 -0.09 0.18 1.001 4199.18 2933.38
is_A:day 0.30 -0.24 0.85 1.002 3657.22 3219.69
is_B:day 0.03 -0.58 0.62 1.002 2777.06 3062.17
is_A:weartime_self_cw NA NA NA NA NA NA
is_B:weartime_self_cw NA NA NA NA NA NA
is_A:barriers_self_cw 0.08 -0.03 0.19 1.000 5788.53 2996.24
is_B:barriers_self_cw 0.01 -0.10 0.11 1.001 4903.70 2641.80
Between-Person Effects
is_A:pushing_self_cb -0.40 -1.15 0.37 1.000 4213.05 3065.63
is_B:pushing_self_cb -0.76 -1.52 0.04 1.001 2795.76 2872.61
is_A:pushing_partner_cb 0.10 -0.63 0.86 1.001 4719.23 2843.64
is_B:pushing_partner_cb 0.51 -0.54 1.59 1.000 2620.35 2881.72
is_A:weartime_self_cb NA NA NA NA NA NA
is_B:weartime_self_cb NA NA NA NA NA NA
is_A:barriers_self_cb 0.10 -0.11 0.32 1.000 3490.33 2884.67
is_B:barriers_self_cb -0.16 -0.46 0.14 1.000 4182.24 3300.63
Random Effects
sd(is_A) 0.16 0.01 0.42 1.00 1611.07 2199.42
sd(is_B) 0.26 0.01 0.64 1.01 980.82 1268.92
sd(is_A:pushing_self_cw) 0.29 0.11 0.48 1.00 1189.28 1251.93
sd(is_B:pushing_self_cw) 0.23 0.04 0.43 1.01 895.24 1039.10
sd(is_A:pushing_partner_cw) 0.13 0.00 0.39 1.00 917.83 1897.62
sd(is_B:pushing_partner_cw) 0.12 0.01 0.35 1.01 1445.05 1996.32
Additional Parameters
ar[1] 0.13 0.00 0.27 1.00 2502.46 2816.40
nu NA NA NA NA NA NA
shape NA NA NA NA NA NA
sderr NA NA NA NA NA NA
sigma 1.00 0.92 1.08 1.00 2258.23 2754.95

Compare Models

loo_compare(loo_ind, loo_apim)
##                        elpd_diff se_diff
## pushing_reactance_ind   0.0       0.0   
## pushing_reactance_apim -6.0       3.5
bayes_R2(pushing_reactance_apim)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.1868663 0.01869625 0.1482543 0.2219012
bayes_R2(pushing_reactance_ind)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.1611657 0.0191676 0.1218981 0.1972937

Report models for pushing side by side

Indistinguishable Models

pushing_all_ind_models <- report_side_by_side(
  pushing_pa_sub_ind,
  pushing_pa_obj_log_ind,
  pushing_aff_ind,
  pushing_anticipAff_ind,
  pushing_reactance_ind,
  
  model_rows_random = model_rows_random_pushing_ind,
  model_rows_fixed = model_rows_fixed_pushing_ind
) 
## [1] "pushing_pa_sub_ind"
## [1] "pushing_pa_obj_log_ind"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pushing_aff_ind"
## [1] "pushing_anticipAff_ind"
## [1] "pushing_reactance_ind"
# pretty printing
pushing_all_ind_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_ind()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pushing_pa_sub_ind 95% CI pushing_pa_sub_ind exp(Est.) pushing_pa_obj_log_ind 95% CI pushing_pa_obj_log_ind b pushing_aff_ind 95% CI pushing_aff_ind b pushing_anticipAff_ind 95% CI pushing_anticipAff_ind b pushing_reactance_ind 95% CI pushing_reactance_ind
Intercept 49.28* [37.64, 64.23] 127.21* [111.99, 144.82] 5.00* [ 4.77, 5.24] 1.91* [ 1.44, 2.37] 0.48* [ 0.25, 0.71]
Within-Person Effects
pushing_self_cw 1.00 [ 0.90, 1.11] 1.00 [ 0.97, 1.04] -0.03 [-0.08, 0.03] 0.01 [-0.11, 0.12] -0.02 [-0.10, 0.07]
pushing_partner_cw 1.02 [ 0.92, 1.13] 0.99 [ 0.96, 1.03] 0.03 [-0.02, 0.09] -0.07 [-0.19, 0.05] 0.02 [-0.07, 0.11]
day 0.86 [ 0.64, 1.15] 1.00 [ 0.88, 1.14] -0.04 [-0.24, 0.15] -0.25 [-0.66, 0.19] 0.27 [-0.12, 0.66]
weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cw 1.02 [ 0.95, 1.08] 0.98 [ 0.96, 1.01] 0.00 [-0.04, 0.03] -0.01 [-0.09, 0.06] 0.05 [-0.03, 0.12]
Between-Person Effects
pushing_self_cb 0.59 [ 0.34, 1.03] 0.82 [ 0.64, 1.07] -0.07 [-0.49, 0.33] -0.34 [-1.22, 0.55] -0.71* [-1.22, -0.20]
pushing_partner_cb 1.50 [ 0.91, 2.49] 0.86 [ 0.66, 1.12] 0.22 [-0.16, 0.60] -0.05 [-0.91, 0.82] 0.27 [-0.31, 0.87]
weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
barriers_self_cb 1.03 [ 0.90, 1.17] 0.99 [ 0.93, 1.05] -0.08 [-0.18, 0.02] 0.17 [-0.06, 0.39] 0.01 [-0.17, 0.18]
Random Effects
sd(Intercept) 0.63 [ 0.44, 0.87] 0.30 [0.22, 0.41] 0.62 [0.47, 0.81] 1.22 [0.92, 1.60] 0.17 [0.01, 0.38]
sd(pushing_self_cw) 0.31 [ 0.13, 0.51] 0.10 [0.04, 0.17] 0.07 [0.01, 0.15] 0.15 [0.01, 0.34] 0.23 [0.11, 0.37]
sd(pushing_partner_cw) 0.28 [ 0.11, 0.48] 0.04 [0.00, 0.11] 0.12 [0.03, 0.22] 0.17 [0.01, 0.39] 0.08 [0.00, 0.23]
Additional Parameters
ar[1] 0.04 [-0.93, 0.95] 0.27 [0.21, 0.33] 0.27 [0.22, 0.33] 0.34 [0.28, 0.39] 0.16 [0.03, 0.28]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.43 [ 0.40, 0.47] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.04 [ 0.00, 0.12] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.53 [0.51, 0.55] 0.89 [0.86, 0.92] 1.88 [1.80, 1.95] 1.01 [0.94, 1.09]
# prepare export to excel
all_models$Pushing_ind <- prepare_for_export(pushing_all_ind_models)

APIMs

pushing_all_apim_models <- report_side_by_side(
  pushing_pa_sub_apim,
  pushing_pa_obj_log_apim,
  pushing_aff_apim,
  pushing_anticipAff_apim,
  pushing_reactance_apim,
  
  model_rows_random = model_rows_random_pushing_apim,
  model_rows_fixed = model_rows_fixed_pushing_apim
) 
## [1] "pushing_pa_sub_apim"
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## [1] "pushing_pa_obj_log_apim"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate = exponentiate, : Coefficients were
## exponentiated. Double check if this was intended.
## [1] "pushing_aff_apim"
## [1] "pushing_anticipAff_apim"
## [1] "pushing_reactance_apim"
# pretty printing
pushing_all_apim_models %>%
  print_df() %>%
  add_header_above(
    c(" ", "Subjective MVPA" = 2, "Device-Based MVPA" = 2, 
      "Mood" = 2, "Anticipatory Affect" = 2, 
      "Reactance" = 2)
  ) %>%
  packing_apim()
Subjective MVPA
Device-Based MVPA
Mood
Anticipatory Affect
Reactance
IRR pushing_pa_sub_apim 95% CI pushing_pa_sub_apim exp(Est.) pushing_pa_obj_log_apim 95% CI pushing_pa_obj_log_apim b pushing_aff_apim 95% CI pushing_aff_apim b pushing_anticipAff_apim 95% CI pushing_anticipAff_apim b pushing_reactance_apim 95% CI pushing_reactance_apim
is_A 64.43* [46.15, 89.96] 129.10* [106.65, 155.78] 5.00* [ 4.71, 5.29] 1.60* [ 0.96, 2.23] 0.30 [-0.01, 0.60]
is_B 35.31* [25.01, 49.86] 126.93* [104.96, 153.25] 4.90* [ 4.59, 5.21] 2.02* [ 1.41, 2.65] 0.74* [ 0.41, 1.08]
Within-Person Effects
is_A:pushing_self_cw 1.05 [ 0.92, 1.21] 0.99 [ 0.95, 1.04] -0.01 [-0.08, 0.07] 0.04 [-0.11, 0.19] 0.00 [-0.12, 0.11]
is_B:pushing_self_cw 0.92 [ 0.79, 1.09] 1.02 [ 0.96, 1.08] -0.06 [-0.14, 0.03] -0.07 [-0.25, 0.10] -0.05 [-0.18, 0.08]
is_A:pushing_partner_cw 1.03 [ 0.89, 1.20] 1.03 [ 0.97, 1.09] 0.00 [-0.08, 0.09] 0.00 [-0.18, 0.17] 0.00 [-0.13, 0.13]
is_B:pushing_partner_cw 1.00 [ 0.87, 1.15] 0.98 [ 0.93, 1.02] 0.05 [-0.02, 0.13] -0.09 [-0.25, 0.06] 0.05 [-0.09, 0.18]
is_A:day 0.73 [ 0.48, 1.12] 1.08 [ 0.91, 1.28] -0.06 [-0.33, 0.21] 0.06 [-0.51, 0.64] 0.30 [-0.24, 0.85]
is_B:day 1.07 [ 0.69, 1.65] 0.94 [ 0.78, 1.12] -0.01 [-0.26, 0.25] -0.55 [-1.13, 0.05] 0.03 [-0.58, 0.62]
is_A:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cw NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cw 1.03 [ 0.93, 1.14] 1.01 [ 0.98, 1.04] -0.02 [-0.07, 0.04] 0.01 [-0.10, 0.12] 0.08 [-0.03, 0.19]
is_B:barriers_self_cw 0.99 [ 0.91, 1.08] 0.98 [ 0.95, 1.00] 0.02 [-0.03, 0.06] -0.03 [-0.13, 0.06] 0.01 [-0.10, 0.11]
Between-Person Effects
is_A:pushing_self_cb 0.47 [ 0.19, 1.15] 0.81 [ 0.56, 1.19] 0.42 [-0.22, 1.08] 1.01 [-0.47, 2.54] -0.40 [-1.15, 0.37]
is_B:pushing_self_cb 0.73 [ 0.35, 1.52] 0.77 [ 0.52, 1.13] -0.28 [-0.84, 0.26] -0.80 [-1.96, 0.41] -0.76 [-1.52, 0.04]
is_A:pushing_partner_cb 1.31 [ 0.68, 2.60] 1.12 [ 0.77, 1.64] -0.07 [-0.61, 0.49] -1.21 [-2.51, 0.05] 0.10 [-0.63, 0.86]
is_B:pushing_partner_cb 2.03 [ 0.84, 4.97] 0.74 [ 0.48, 1.16] 0.41 [-0.23, 1.06] 1.16 [-0.36, 2.58] 0.51 [-0.54, 1.59]
is_A:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_B:weartime_self_cb NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA
is_A:barriers_self_cb 1.08 [ 0.89, 1.31] 1.07 [ 0.99, 1.17] 0.01 [-0.14, 0.16] 0.21 [-0.11, 0.54] 0.10 [-0.11, 0.32]
is_B:barriers_self_cb 1.04 [ 0.82, 1.30] 0.98 [ 0.89, 1.09] -0.07 [-0.25, 0.11] 0.04 [-0.35, 0.42] -0.16 [-0.46, 0.14]
Random Effects
sd(is_A) 0.63 [ 0.39, 0.97] 0.34 [0.24, 0.47] 0.79 [0.61, 1.02] 1.43 [1.08, 1.88] 0.16 [0.01, 0.42]
sd(is_B) 0.47 [ 0.29, 0.68] 0.34 [0.23, 0.48] 0.68 [0.51, 0.89] 1.52 [1.15, 1.98] 0.26 [0.01, 0.64]
sd(is_A:pushing_self_cw) 0.12 [ 0.01, 0.33] 0.09 [0.01, 0.18] 0.07 [0.00, 0.18] 0.15 [0.01, 0.37] 0.29 [0.11, 0.48]
sd(is_B:pushing_self_cw) 0.46 [ 0.14, 0.85] 0.11 [0.01, 0.20] 0.12 [0.01, 0.27] 0.19 [0.01, 0.49] 0.23 [0.04, 0.43]
sd(is_A:pushing_partner_cw) 0.43 [ 0.11, 0.84] 0.11 [0.02, 0.20] 0.17 [0.04, 0.32] 0.28 [0.02, 0.60] 0.13 [0.00, 0.39]
sd(is_B:pushing_partner_cw) 0.11 [ 0.01, 0.30] 0.03 [0.00, 0.09] 0.10 [0.01, 0.24] 0.11 [0.00, 0.31] 0.12 [0.01, 0.35]
Additional Parameters
ar[1] 0.04 [-0.92, 0.94] 0.21 [0.14, 0.28] 0.18 [0.12, 0.24] 0.24 [0.18, 0.30] 0.13 [0.00, 0.27]
nu NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
shape 0.43 [ 0.40, 0.47] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sderr 0.05 [ 0.00, 0.13] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA] NA [ NA, NA]
sigma NA [ NA, NA] 0.52 [0.50, 0.54] 0.85 [0.82, 0.89] 1.80 [1.73, 1.88] 1.00 [0.92, 1.08]
# prepare export to excel
all_models$Pushing_apim <- prepare_for_export(pushing_all_apim_models)
openxlsx::write.xlsx(all_models, file = "Output/AsPreregistred/AllModels.xlsx")